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I.  INTRODUCTION 

The  calculation  of  viscous  flow  fields  past  high  speed  axi- 
symmetric  blunt  bodies  is  of  prime  interest  to  the  designer  of 
the  reentry  space  vehicles.  The  variety  of  conditions  met  during 
a vehicle  reentry  requires  the  solution  to  be  valid  over  a wide 
range  of  Reynolds  number,  ranging  ‘from  a low  Reynolds  number  at 
high  altitude  to  a high  Reynolds  number  at  low  altitude. 

p 

For  high  Reynolds  number  flows  past  blunt  bodies,  a number  of 
effective  methods  have  been  developed  for  analyzing  the  inviscid 
flow  and  the  associated  boundary  layer  flow  over  the  body.  In 
flow  regimes  of  intermediate  to  low  Reynolds  numbers,  the  viscous 
region  incompasses  a significant  portion  of  the  shock  layer  and  the 
classical  boundary-layer  approaches  can  not  be  used.  The  idea 
of  using  the  second  order  boundary  layer  equations  [1]  to 
compute  the  flow  field  becomes  appealing  in  such  cases;  however, 
this  approach  leads  to  several  computational  difficulties. 

Aside  from  the  excessive  computing  time,  one  experiences  problems 
with  the  strong  vorticity  interaction  region  which  causes 
difficulty  in  the  matching  procedure  between  the  viscid  and 
inviscid  region. 

Because  of  the  difficulties  mentioned  above,  it  is  desirable 
to  seek  an  alternate  method  of  approach  to  the  problem.  The 
use  of  the  full  Navier-Stokes  equations  [2,  31  has  been  successful 
in  providing  solutions  for  the  stagnation  regions  but  generally 
have  been  applied  for  only  one  to  two  nose  radius  downstream. 

The  complexity  of  the  solution  procedures  of  the  Navier-Stokes 
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equations  restrict  their  applications  especially  in  the  downstream 
direction.  Because  of  this  difficulty,  attention  has  lately 
been  turned  toward  the  viscous  shock  layer  techniques.  The 
governing  equations  of  these  techniques  are  basically  parabolic 
in  nature  and  their  numerical  solution  is  a straightforward 
streamwise  marching  procedure. 

The  viscous  shock  layer  equations,  which  were  developed 
primarily  by  Davis  [4] , contain  all  of  the  terms  in  the  Navier 
Stokes  equations  which  contribute  to  second  order  effects  in 
both  the  inviscid  and  viscid  regions.  If  the  pressure  gradient 
normal  to  the  body  surface  is  assumed  to  be  established  entirely 
by  centrifugal  effects,  the  thin  shock  layer  version  of  the 
equations  is  obtained.  While  the  solution  of  the  thin  layer 
equation  is  straightforward  with  techniques  similar  to  the 
boundary  layer  marching  techniques,  the  solution  of  the  full 
shock  layer  has  encountered  difficulties.  In  an  attempt  to  solve 
the  full  layer,  Davis  developed  a method  wherein  the  solution 
is  obtained  through  a relaxation  process  from  the  thin  shock 
layer  solution.  While  this  method  was  successful,  it  is  sensitive 
and  can  encounter  divergent  behavior  in  the  relaxation  scheme 
used  to  relieve  the  thin  shock  layer  assumptions  [5].  This 
difficulty  is  most  severe  for  slender  bodies.  In  an  attempt 
to  overcame  this  problem,  a different  relaxation  technique  was 
developed  by  Werle  et  al.  [5].  This  technique  utilized  an 
artificial  time  coordinate  to  relax  the  shock  from  an  initial 
shape.  Limited  success  has  been  reported  when  under-relaxation 
factors  were  used  to  remove  instabilities  that  occurred  in  the 


6 


AEDC-TR-79-25 


solution  far  downstream  from  the  stagnation  point. 

The  Werle  et  al.  [5]  and  Davis  {4]  methods  are  both  essen- 
tially relaxing  the  shock  shape  to  take  into  account  downstream 
influences  and  consequently  acknowledge  the  boundary  value 
nature  of  the  problem.  In  both  cases,  difficulties  were 
encountered  whenever  the  shock  layer  thickness  is  large  and 
especially  far  downstream  on  blunt  slender  bodies  where  the 
inviscid  region  encompasses  a significant  portion  of  the  total 
shock  layer  thickness.  These  difficulties  do  not  stem  from  the 
shock  shape  relaxing  technique  but  rather  occur  due  to  other 
fundamental  reasons. 

It  is  known  from  small  disturbance  theory  in  hypersonic 
flow  [6]  that  the  streamwise  velocity  variation  does  not  affect 
the  solution  of  the  other  flow  variables  since  the  group  of 
equations  representing  the  continuity,  normal  momentum  and 
energy  are  completely  uncoupled  from  the  longitudinal  momentum. 
Therefore  a similar  behavior  of  the  shock  layer  equations  is 
expected  in  the  inviscid  region  far  downstream  on  a slender 
blunt  body. 

Far  downstream  on  a slender  body,  the  boundary  layer  is 
comparatively  thin  and  the  flow  in  the  shock  layer  is  pre- 
dominately inviscid.  A general  description  of  the  inviscid 
flow  could  be  obtained  from  the  solution  of  the  inviscid  shock 
layer  equations.  For  slender  bodies  with  relatively  small 
surface  slope,  the  approximations  of  hypersonic  small  distur- 
bance theory  is  quite  important.  If  these  approximations  are 
applied  to  the  inviscid  shock  layer  equations,  and  only  first 
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order  terms  are  kept,  the  normal  momentum,  continuity  and  energy 
equations  become  uncoupled  from  the  longitudinal  momentum 
equation.  It  is  therefore  clear  that  in  this  case  the  solution 
of  the  normal  momentum,  continuity  and  energy  equations  will 
not  depend  strongly  on  the  solution  of  the  tangential  momentum 
equation . 

All  numerical  methods  for  handling  the  shock  layer  equations 
solve  the  governing  equations  in  a successive  manner  where  the 
solution  of  the  tangential  momentum  equation  drives  the  solution 
of  the  normal  momentum  and  continuity  equations.  These  numerical 
methods  become  improper  for  the  flow  far  downstream,  since  the 
tangential  momentum  equation  does  not  have  a predominant 
influence  on  the  solution.  A more  adequate  solution  technique 
is  to  solve  the  equations  simultaneously;  this  will  ensure  the 
proper  coupling  between  each  equation  in  the  development  of  the 
solution . 

In  reference  [5],  difficulties  in  the  solution  were  reported 
as  occurring  in  the  iterative  loop  that  is  used  to  solve  the 
continuity  and  the  normal  momentum  equations.  This  difficulty  was 
reported  as  an  oscillatory  behavior  of  the  solution  from  one 
iteration  to  next.  It  was  suggested  in  the  same  reference 
that  the  continuity  and  normal  momentum  equations  should  be 
solved  simultaneously  as  a coupled  set  to  overcome  this  mis- 
behavior. Waskiewicz  et  al.  [7]  have  applied  this  suggestion  and 
have  reported  good  improvement  in  the  solution  obtained  on  a 
sphere.  It  seems  that  the  coupling  has  strengthened  the  mutual 
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dependence  of  the  continuity  and  normal  momentum  equations. 

In  general,  better  solutions  are  expected  if  all  the  govern- 
ing equations  are  solved  simultaneously  as  a coupled  set,  as 
the  interaction  and  dependence  of  the  equations  on  each  other 
will  be  automatically  guaranteed  and  the  physical  nature  of  the 
problem  will  be  automatically  satisfied.  While  solving  the 
equations  as  a coupled  set  has  shown  excellent  results  in 
handling  inviscid  shock  layer  equations  [8,  9],  it  has  not  yet 
been  used  in  the  viscous  shock  layer  equations . The  aim  of 
the  present  study  is  to  develop  a numerical  algorithm  that  is 
capable  of  solving  two  first  order  and  two  second  order 
equations  simultaneously  and  to  assess  its  applicability  to  the 
solution  of  the  full  viscous  shock  layer  equations. 
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II.  GOVERNING  EQUATIONS 

The  viscous  shock  layer  equations  used  in  this  study  have 
been  developed  by  Davis  [4].  A detailed  derivation  of  these 
equations  are  given  in  references  [4]  and  [10],  and  they  are 
only  summarized  here.  The  Navier  Stokes  equations  are  first 
written  in  a boundary  layer  coordinate  system  (see  Fig.  1) . 

Terms  higher  than  second  order,  in  the  inverse  square  root  of 
the  Reynolds  number,  from  both  a viscid  and  an  inviscid  view- 
point are  then  neglected.  The  resulting  equations  are: 

Continuity : 

[ (r+ncos4>)  ^ pu]  + ( (l+<n)  (r+ncos4 ) pv]  = 0 (la) 

Longitudinal  Momentum: 

p{u  u /(l+Kn)  + v u + kuv/ ( 1+icn)  } + p /(l+ien) 
s n s 

= [e2/  (1+tcn)  2 (r+ncos40  ] [ (1+Kn)  2 (r+ncos<f>)  ^t]  (lb) 

n 

where, 

t = u [u  - tcu/(l+Kn)]  (1c) 

n 

Normal  Momentum: 

p{u  v /(1+icn)  + v v - ku2/{1+  n)  } + p =0  (Id) 

s n n 

which  with  the  thin  shock  layer  approximation  becomes, 

- p ku2/( 1+Kn)  + pn  = 0 (le) 
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Energy  Equation: 

2 2 

p{uT  /(l+Kn)  + vT  } - u p /(l+<n)  - v p = e t /y 
s n rs  rn 


+■  [e2/ (l+<n)  (r+ncos<jj) -1  ] [ (1+icn)  {r+ncosip)  ^q]n  (If) 

where 

q = y T /a  (lg) 

Equation  of  State: 

p = ( Y- 1 ) PT/Y  (lh) 

Viscosity  Law: 

u = T3/2  (H-c')/(T+c')  (li) 


where 

c’  = c*/M2  T*  (y-1) 

* 

and  c is  taken  to  be  198. 6°R  for  air. 

The  surface  boundary  conditions  are  the  no  slip 


conditions : 

u(s,o)  = v(s,o)  = 0 

(2a) 

and 

T(s,o)  = T 

w 

(2b) 

while  at  the  shock  location  the  oblique  shock  relations  are  used. 
These  relations  are  given  as: 

u . = u . sin(a+6)  + v , cos(a+£)  (3a) 

sh  sh  sh 

vgh  = - ugh  cos(a+S)  + vgh  sin (a+B)  (3b) 
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where  u . and  v , are  velocity  components  at  the  shock  given 
sh  sh 

as 


cosa 

- sina/p 


sh 


and 


psh 

Psh 


YPsh/(V-1)Tsh 

(2/(y+1)  1 sin^a  - (y-1)  /Y  (Y+D  M2 
(ush  * cosa)2/2  + { [4Y/(Y+D2Isin2a  +[2/(y-1) 
- 4(y-1)/(Y+1)2]/M2  - 4/{y+1)2  M-i  sin2a>/2 


(3c) 

(3d) 


( 3e) 
( 3 f ) 


(3g) 


For  ease  in  the  numerical  computations  and  also  for  shock 
relaxing  considerations,  the  above  equations  are  normalized  by 
using  the  following  variables: 


" = n/nsh 

K - s 

5 * u/ush 

? - v/vsh 

■ 

t = T/Tsh 

p = p/psh 

5 - o/Psh 

“ - “■'“ah 

(4a-h) 

The  resulting  equations  are  two  second  order  and  two  first 
order  differential  equations  and  are  given  by: 

s-Momentum  Equation: 

32u/3n2  + c^Ou/an)  + a2u  + a3  + a4(3u/30  = 0 (5) 


where 
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p , u , n , n , --  P u v u n v.  r- 

sh  sh  sh  sh  pun  _ sh  sh  sh  pv 

2 1+Kn  . n - 2 

E v,  sh  U e 


<n  . cosebn  , 

+ y /V  + t + sh 

n ■L+Knshn  Y+nshncos^ 


a2 


p , u , n , n , — o . v ' n . <n  . 

sh  sh  sh  sh  pu  _ sh  sh  sh  sh  - _v_ 

2 1+Kn^.  n - 2 l+<n  , n 

e ySh  Sh  y E Ush  sh  U 


n , icn  , cos$n  , <n  , 

- K 5h  u /u  Sh  + Sh  ) y ( Sh 

1+Kn  . n it  '1+Kn  , n r+n  , ncosi  ‘ l+<n  , n 
sh  sh  sh  sh 


a3  = - 72- 


p , n , n , 
sh  sh  sh 


n u n p , 

, ■■sh  1/y  1/u  . (pr  - r, p + -^h  p} 
l+<na.  sh  neV,  n PcV,  w 


a4  = - (Bshushnsh/e2u5h)  (nsh/tl+,<nshr!)  “ 


Energy  Equation: 


32t/3n2  + ax  Ot/3n)  + fi2  t + a3  + a4  { 3 t/9C ) " 0 


where 


P u u , n a n , --  o , v , n , a 

r _ sh  sh  sh  sh  oun  _ sh  sh  sh  ,,  pv 

1 2 1+Kn  . n - ” 2 - 

e ush  sh  V € ush  y 

ten  . cosdrn  , 

+ y /u  + , Sh  ■ + 

n 1+Knghn  r+nshncos4 


“2  = -(BshushTih/B2"shTsh,<nsh/(1+<n-r’,)  * ^ 
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P u“  n , U 

sh  sh  sh  [„_sh__ 


£ >‘.hT»h 


^+KnshT1  n 


SK  ^ 


sh 


2 

v u . a _ 

+ vS]  + ( u - 

u . rn  n 


Kn 


sh 


sh 


sh 


1+Knshn 


u>2 


(6c) 


54  = - <°^h^hWe\h),nsh/tl+'(nsht',)  f 


( 6d) 


Continuity  Equation: 


Insh(r+nshnoos,t)  Sshush  5S|5  + l(rtVcos*|1‘ 


U1+™  n)pghvsh  p?  - nsh  cshush  pun } ] „ » 0 


(7) 


n-Momentum  Equation: 


TTS^nT  <’c  - »>*  -n  + ^ + ^ ^ % 


PU 


u 


1+*nsh  vsh 


sh  --2 

pu  + 


sh 


f)shushVshnsh 


P„  - o 


(8) 


where  with  the  thin  layer  approximation  this  equation  becomes, 


p = [u/(l«nshn>]<Pshu^nsh/Pah>  in2 


(8a) 


The  remaining  equations  are: 
The  Equation  of  State: 

P = Pt 


(9) 
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and  the  viscosity  law, 

V - [(Tgh  + c’)/(Tsht  + c')]  t3/2 


(10) 


The  surface  boundary  conditions  are 

u = 0,  v = 0,  and  t = t (11) 

and  the  conditions  at  the  shock  ares 

u = v = t = p=  p = y = l (12) 


An  equation  of  global  mass  conservation  can  be  obtained  from 
equation  (7)  by  integrating  from  r = 0 to  n = 1 while  holding 
5 constant.  This  results  in 


dm 

d5 


(r+n  , cos4) 
sh 


[nshpshush 


<1+Knsh)pshvsh1 


(13a) 


where 


1 

m = / nsh(r+nshr,COS?)  °shush  P"  dn  (13b) 

is  proportional  to  the  rate  of  mass  flux  between  the  body  and 
shock  at  a given  position  on  the  body  surface.  A limiting  form 
of  equations  (13a)  and  (13b)  is  obtained  at  the  stagnation  point 
by  applying  the  condition  5=0  and  using  L'Hopitals  Rule. 

This  limiting  form  is  used  to  determine  the  shock  stand-off 
distance  at  the  stagnation  point. 

Equations  (5)  to  (10)  constitute  the  complete  set  of  governing 
equations  for  the  unknowns  u,  v,  t,  p,  p and  p.  These  equations 
are  solved  along  with  the  surface  conditions  given  by  equation  (11) 
and  the  shock  conditions  given  by  equation  (12) . 
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III.  THE  "UNSTEADY " SHOCK  LAYER  EQUATIONS 

Werle  et  al.  [5]  have  developed  a time  relaxation  scheme 
where  an  initial  guess  of  the  shock  shape  is  relaxed  in  an 
artificial  time  like  manner  toward  the  "steady  state"  solution. 

The  scheme  was  developed  to  overcome  the  divergent  behavior  in 
the  iteration  scheme  developed  by  Davis  [4],  Two  versions,  but 
essentially  the  same,  have  been  developed  in  which  the  first 
relaxes  the  shock  wave  thickness  [5]  , while  the  second  relaxes 
the  shock  radius  [11].  The  first  approach  is  very  suitable  and 
directly  handles  the  variable  which  appears  always  in  the 
equations.  But  the  second  is  more  general  and  it  could  handle 
cases  that  involve  discontinuity  in  the  surface  curvature. 

For  reasons  of  generality,  the  second  version  will  be  adopted 
and  will  be  summarized  here. 

In  order  to  demonstrate  the  adopted  time  relaxation  technique, 
the  "unsteady"  shock  layer  equations  should  first  be  obtained. 

In  reference  [5],  Werle  et  al.  have  reasoned  that  an  artificial 
time  term  can  be  added  to  the  equations  through  the  shock 
variable  derivatives.  For  this  reason,  the  governing  equations 
have  to  be  in  a normalized  form  to  provide  explicitly  the  shock 
derivatives.  In  Appendix  A,  the  shock  derivatives  are  -given  in 
terms  of  the  shock  slope  (da/ds) . Since  da/ds  represents  the 
curvature  of  the  shock  shape,  a relation  between  da/ds  and 
the  second  derivative  of  the  shock  radius  is  derived  in 
Appendix  B and  is  given  by 

< 2 

da  _ d2R  , cos2  ( ct— 4> ) . dR  ,icsin  ( 2a-2a)  + nsh<  cos  ^ , 
ds  ds2  l(l+<nsh)cos^  " ds  coscKl+Kn^)  cos<>  (1+7nsh?_  (14) 
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This  relation  differs  from  that  of  reference  [11]  in 

including  the  effect  of  the  variation  of  the  body  curvature  (<'). 

The  artificial  time  like  term  is  then  introduced  to  give 

» 2 , 

da  _ .32R  9R,  . cos2(a-j>)  dR  .xsin  ( 2ot—  24> ) + nshIC  COS  a . 

ds  gg2  " 3tJ  l(l+Kns)cos<},  “ ds  lcos*(l+Kns)  COS!j)  d+Knsh)  2 

(15) 


At  the  stagnation  point  <}>  = 90°  and  this  equation  becomes 
singular.  However,  if  the  limit  is  taken  of  this  equation 
as  s tends  to  zero,  the  following  equation  is  obtained 


32n 


[- 


sh 


da  1 

E = ^ ^ 


1 !^sh 

2 3t 


] - 1 


(16) 


Upon  substitution  for  da/ds  in  the  shock  derivative  relations 
given  in  Appendix  A and  using  them  in  the  governing  equations, 
the  "unsteady"  equations  are  achieved. 

I 

In  summary,  equations  (5)  to  (8)  together  with  equations 
(Al,  A2 , A3,  A4 ) and  equation  (15)  represents  the  set  of 
11  unsteady"  shock  layer  equations.  The  set  is  a form  suitable  for 
numerical  computations.  However,  for  time  relaxation  considerations, 
it  is  necessary  to  write  the  longitudinal  momentum  in  its  final 
form  given  by 


32u 


H + *1  iw  + M-H  - ffl  + 63  ft  + *4  - 85  ~ = o 

d T]  o S 


• 3 R 3Rn 


3R 


3u 


(17) 


where  8^  B2,  3-j,  S4  and  35  are  obtained  from  reference  [11] and 
are  listed  here  in  Appendix  C. 
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IV.  NUMERICAL  ANALYSIS 


1.  General  Considerations; 

The  "unsteady"  governing  equations  are  solved  using  a time 
relaxing  approach  equivalent  to  the  implicit  alternating  direction 
method.  This  approach  utilizes  two  steps  in  which  the  first 
yields  the  flow  properties  in  the  shock  layer  while  the  second 
step  is  to  update  the  shock  shape  itslef.  This  is  demonstrated 
by  writing  the  "unsteady"  longitudinal  momentum  equation  (17) 
in  two  time  steps: 


First  Step  (Star  Sweep) 


_ .n  * .n  , At 
from  t to  t = t + ~— 


32u*  * 3u*  * > 32Rn 

7 


3R*.  , p*  5Rn  , , p* 

3t”]  + S3  + B4  + B5  35 

At 


= 0 (18) 


* n+1  * 

Second  Step  (Final  Sweep) : from  t to  t = t + 


3R 


* * 

62  -TT-  - s2  at 

as 


n+1 


+ 8 


* 3R^ 
3 3s 


n+1 


a2- 

[2-S  + 

3n 


S — 
5l  3n 


+ 65  If  + 64'‘  * 


(19) 


Using  equation  (C6)  of  Appendix  C,  equation  (19)  can  be 
written  in  a form  independent  of  n by  using  equation  (18)  to 
give 


32Rn+1 

3s2 


n a_n+l 

[2k  tan  (ct-4>)  + 3i 

sn 


— nn+1 
At  R 


..  2 _n  n , aRr 

H-  + [2K  +ifnrr!  sir 

3 s sh 


+ 


2_ 

At 


* 

(2R 


Rn) 


(20) 
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The  boundary  conditions  associated  with  the  star  sweep  of  the 
unsteady  equations  are  typical  no  slip  conditions  (11)  at  the 
surface  and  Rankine— Hugoniot  conditions  at  the  shock  location  (12)  . 
The  boundary  conditions  associated  with  the  final  sweep  are  the 
same  as  those  used  in  reference  [11]  and  are  given  as: 


at  s 
at  s 


0 

s 

max 


,n+l 


^n+1 

max 


0 

* 

R 

max 


(21a) 

(21b) 


2 • Method  of  Solution  of  the  Star  Sweep: 

The  deviation  point  between  the  present  work  and  that  of 
references  [4]  and  [5]  is  the  way  the  equations  are  solved  during 
the  star  sweep.  Previously,  the  equations  were  solved  in  a 
successive  manner  where  the  solution  of  one  equation  drives 
the  solution  of  the  other.  This  technique  is  widely  used  in 
many  numerical  calculations  and  is  known  as  the  cascading  scheme. 
In  the  present  work,  the  equations  will  be  solved  simultaneously. 
This  technique  will  be  referred  as  the  coupling  scheme.  In 
reference  [7]  a combination  of  the  two  schemes  were  used  where 
the  two  second  order  equations,  s-momentum  and  energy,  are  solved 
in  a successive  manner  while  the  continuity  and  n-momentum  are 
solved  simultaneously.  This  combination  has  indicated  an  improve- 
ment in  the  solution  as  discussed  earlier  in  the  report. 

The  equations  to  be  solved,  during  the  star  sweep,  are  the 
s-momentum,  energy,  n-momentum  and  the  continuity.  When  the 
equation  of  state  is  used  to  replace  the  density  variable  by  the 
pressure  and  temperature  variables,  the  governing  equations  are 
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mainly  two  first  order  and  two  second  order  equations  in  the  un- 
knowns u,  v,  p and  t.  They  require  six  boundary  conditions,  three 
on  the  body  surface  and  three  at  the  shock.  Actually,  there  are 
seven  boundary  conditions  available,  four  at  the  shock  and  three 
at  the  body  surface,  but  there  is  an  additional  unknown  n^-  Thus 

the  equations  and  the  boundary  conditions  form  a complete  system  of 

* 

equations.  In  this  solution  procedure,  nsh  is  assumed  to  be  known 
and  is  iterated  on  to  satisfy  the  extra  boundary  condition.  In 
this  case,  the  extra  boundary  condition  is  chosen  to  be  v(n=l)  = 1. 

v 

The  governing  equations  (5) -(8)  are  nonlinear  equations. 

It  is  first  convenient  to  linearize  the  governing  equations  with 
the  following  relations,  where  the  quantities  with  superscript  ° 
are  evaluated  from  previous  iteration 


XY 

=5 

X«Y  + Y° X - X°Y° 

(22a) 

XYZ 

= 

Xe>YaZ  + Y°Z°X  + Z°X°Y  - 2 

X°Y°X° 

(22b) 

XZ 

Y 

- 

Xa  . 2°  X°  Z° 

F Z Y“  X Y«2  Y 

(22c) 

Where  X, 

Y 

and  Z can  represent  any  of 

the  variables  u. 

v,  p and  t 

and  their  derivatives.  If  X° , Y°,  Z°  values  are  taken  from 
previous  £ solution,  the  linearization  relations  are  second  order 
in  AS  and  nonlinearity  iterations  are  not  necessary  (This  method 
is  known  as  quasi-linearization) . 

The  resulting  linear  equations,  after  dropping  the  super- 
script , are: 

S-Momentum: 

.2 

+ (a^u^  + a3p^}  + (b^  + b^) 

3 n 

+ (c^u  + C2 v + c^p  + c^t)  + d = 0 (23) 
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Energy : 

a2t 


3n 


2 + U3PC  + a4fc^  + (biun  + b3Pn  + b4V 


+ (c^u  + c^v  + c^p  + c^t)  + d = 0 


(24) 


n- Moment  tain : 

a2V£  + ^b2Vn  + b3Pn*  + ^C1U  + C2V  + C3P  + C4b^  + d = 0 (25> 
Continuity: 

(a  + a3pf  + a4fc^'  + (biun  + b2Vn  + b3pn  + b4tr) 

+ (c^u  + c^v  + c^p  + c^t)  + d = 0 (26) 


The  coefficients  a^,  a2,  ...  are  evaluated  independently  for  each 
individual  equation  and  are  given  in  Appendix  D.  To  solve  these 
equations  numerically,  it  is  necessary  to  write  these  equations 
in  a finite  difference  form.  Let  the  subscript  m denote  the 
station  measured  along  the  body,  while  the  subscript  n denotes 
the  station  measured  normal  to  the  body,  the  derivatives  in  the 
r direction  are: 


<£> 


m,  n 


= [An  ./An  (An  + An  .)]  w , 
n-l  n n n-l  m,n+i 


+ [<AT'„-iT>n-l)/AVl1n-l1  Wm, n [ A VAnn-l { Ann  + Atln-i  > 1 


W “ ~ An  Ar  1 

m,n-l  6 m>n  n n-l 


(27) 
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= [2/An  (An  + An  , ) 3 v,  - (2/An  An  ,)  w 

„ 2 _ _ n n n-l  in,n+l  n n-1  m,n 

d n ni , n 

1 a4 

+ [2/An  n (An  +An  ,)w  , - ( — x)  An  An  i 

' n-l  n n-l  m,n-l  12  \ 4'  n n-l 

o n m , n 

3 4 

- i (^~)  (An_-Ann  ,)  (An  —An  ,)2 

3 an3  m,n  n n_1  12  3n4  m,n  n n_1  (26) 


While  the  derivative  in  £ direction  is  a two  point  difference 


given  by: 


— = (w  „ - wm  . )/A£ 

o Cf  in  i n m-  / n 


The  n derivatives  are  of  second  order  accurate  in  the  step 
size  if  constant  An  is  chosen.  If  the  £ derivative  is  evaluated 
at  mid  point  (m-^n)  and  other  terms  are  averaged,  one  obtains 
a Crank-Nicolson  scheme.  If  a backward  difference  is  used 
for  9w/3£  at  the  (m,n)  a purely  implicit  scheme  is  generated. 

The  previous  difference  quotients  are  suitable  for  the 
second  order  s-momentum  and  energy  equations.  For  the  first  order 
n-momentum  and  continuity  equations,  it  is  convenient  to  use 
two  point  difference  relations  for  both  n and  £ derivatives  as 


follows 


, “ <wm,n  - “n,n-l)/4Vl  + 0(4,,n-l) 

iuf  n- 


V W f 

Tr  = (W_ 


. - w . , )/A£ 

jn-1!  m-ljn-ij 
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If  the  £ derivative  is  evaluated  at  (m-^n-^)  and  •*—  is  averaged 

dri 

between  the  stations  m-1  and  m,  the  box  scheme  is  obtained. 

If  the  linearized  forms  of  s-momentum  (23),  energy  (24) , 

n-momentum  (25}  and  continuity  (26)  are  evaluated  at  (m-4,n), 

► 

(m-%,n)  , (m-^n-'i)  and  (m-4,n+^)  respectively,  the  resulting 

difference  equations  will  be  in  a form  suitable  for  solution 
with  an  algorithm  similar  to  Thomas  algorithm.  Using  the 
difference  quotients  {27)  — (29)  and  (30)-(31),  these  difference 
equations  are: 


(allu„-l  + + <bllan  + b12vn  + b13pn  + lu‘n> 


* (cUun+l  + c13pn+l)  ” dl 


where  n = 2,  3,  ....  (N-l) . 

The  coefficients  in  the  above  equations  vary  from  one  grid 
point  to  another  and  are  given  in  Appendix  E,  These  linear 
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algebraic  equations  are  to  be  solved  with  the  following  boundary 
conditions , 


u 


1 


0 


v 


1 


0 


t 


1 


w 


“ 1 


t 


N 


1 


(36a) 
(36b) 
(36c) 
(37a) 
(37b) 
(37c) 
( 37d) 


The  difference  equations  together  with  the  boundary  conditions 
lend  themselves  to  solution  by  the  following  recursion  relations: 

u,=D  U+E  v+F  p + G,  t+H  (38a) 

n+1  un+l  n un+l  n Un+X  n un+l  n Un+1 

v ,,  = D u+E  v+F  p + G„  t + H (38b) 

n+1  Vn+1  n vn+l  n vn+l  n Vn+1  n Vn+1 

P+=D  u+E  v+F  p+G  t+H  nn-, 

n+1  pn+l  n pn+l  n pn+l  n pn+l  n Pn+1  (38c) 


t . , = D.  u + E v+F,  P + G.  t+H  (38d) 

n+1  ^+1  n fcn+l  n Vl  n fcn+l  n fcn+l 


The  two  subscripts  marking  the  coefficients  D,  E,  F,  G and  H 
denote  to  which  variable  and  grid  point  they  belong.  The  solution 
steps  using  these  recursion  formulas  are  given  in  Appendix  F . 
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3 . Solution  of  the  Final  Sweep  Equation; 

The  governing  equation  of  this  step,  at  a fixed  time,  is  of 
the  form 


where 


d2Rn+1  A dRn+1 
— 5! + a 


dS  ‘ 


1 dS 


+ a2  Rn+1  + a3  = 0 


(39) 


al  = “[2k  tan(a“*)+  jg— ] 

sn 


a2  = - 2/At 
92Rn 


"sh*' 


3Rn  , 2 


“3  = - +[2«  + iSnT’  + ft  (2R  - 

SS  sn 


with  the  boundary  conditions 
Rn+1 {0 ) » 0 


Rn+^  fSmax)  = R*(Snax) 


The  problem  is  a two  point  boundary  value  problem  and  its 
numerical  solution  is  straightforward  since  it  reduces  to  a 
tridiagonal  form  once  the  following  second  order  difference  quotients 
are  substituted  into  the  differential  equation  (39): 


(Vl  - 2Rm  + VlJ/  aS‘ 


(40a) 


dR 

dS 


(R. 


m+1 


R ,)/2AS 
m- 1 


(40b) 
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4 . Overall  Method  o£  Solution: 

The  overall  method  of  solution  employed  is  as  follows.  An 

initial  guess  is  first  made  for  the  shock  shape  over  the  region 

of  interest.  This  guess  can  be  made  arbitrary  and  the  most 

simple  one  is  to  assume  the  shock  lies  parallel  to  the  body 

surface  with  a constant  shock  thickness  . Based  on  this  guess, 

the  first  and  second  derivatives  of  the  shock  radius  are 

evaluated.  The  star  sweep  equations  are  then  solved  by  starting 

at  the  stagnation  point  where  the  governing  equations  are 

reduced  to  ordinary  differential  equations.  With  initial  guesses 

* 

for  the  flow  profiles  and  for  nsh,  all  the  linearized  governing 

equations  are  solved  simultaneously  with  the  numerical  algorithm 

* 

presented  earlier  to  obtain  the  new  flow  profiles  and  n 

n , at  each  iteration  is  calculated  from  the  equation  repre- 
ss 

senting  the  global  mass  conservation.  The  coefficients  of  the 

linearized  equations  are  then  reevaluated.  Repetition  of  the 

above  steps  is  continued  until  a converged  solution  is  obtained 

at  the  stagnation  point.  The  method  then  steps  along  the  body 

surface.  The  previous  station  values  of  the  flow  profiles  are 

used  at  each  new  step  to  evaluate  the  coefficients  of  the  quasi- 

linearized  equations.  Their  solution  is  then  obtained  by  the  same 

numerical  algorithm.  Iterations  on  the  nonlinear  terms  is  not 

necessary  since  the  solution  will  be  at  least  of  second  order 

accurate  in  A£.  However,  local  iterations  at  each  £ station 

* 

are  still  necessary  to  obtain  the  proper  nsh  which  would  satisfy 
the  extra  boundary  condition  v(ri=l)  = 1 (or  vtn^l)  = vg^)  . 
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A Newton- Raphson  iteration  procedure  is  utilized  to  obtain. the 

'* 

shock  distance  n 

sh 

Once  the  above  procedure  has  marched  over  the  entire  mesh 

n+X  , 

the  final  sweep  is  then  solved  for  ns^  . No  iteration  of  the 
final  sweep  equation  is  required  since  the  equation  is  linear. 
The  shock  shape  solution  obtained  is  then  used  to  initialize  the 
following  star  sweep  in  time.  This  process  is  continued  until 
two  alternate  final  sweeps  converged  to  a desired  degree  of 
accuracy.  The  solution  obtained  is  the  required  "steady  state" 
solution. 
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V.  RESULTS  AND  DISCUSSION 


The  numerical  algorithm  developed  in  Section  IV  was  applied 
to  obtain  the  solution  of  the  full  viscous  shock  layer  equations 
over  blunt  body  configurations.  The  success  of  the  algorithm 
lies  in  its  ability  to  extend  the  flow  field  calculations  far  aft 
of  the  nose  region  without  the  oscillatory  behavior  reported 
in  references  [5]  and  [11].  Numerical  solutions  were  successfully 
generated  up  to  twenty  times  the  nose  radius  downstream  over 
hyperbolic  blunt  bodies.  The  algorithm  was  also  capable  of 
handling  flow  fields  around  slender  bodies  and  this  was  demon- 
strated by  obtaining  solutions  over  hyperboloid  bodies  with 
small  asymptotic  angles.  The  reliability  of  the  method  was 
checked  through  comparisons  with  other  numerical  calculations 
and  with  available  experimental  data. 

In  all  cases,  the  gas  was  assumed  to  be  a perfect  gas  with 
a constant  specific  heat,  constant  Prandtl  number  and  the 
viscosity  was  assumed  to  be  given  by  the  Sutherland  law.  The 
first  case  considered  was  that  of  flow  over  a 22.5°  half  angle 
hyperboloid.  The  flow  conditions  were  the  same  as  those  con- 
sidered in  reference  [5]  with  a free  stream  Mach  number  of 
= 21.75,  free  stream  temperature  of  Tm  _ 351>8oH  wall  to 

stagnation  temperature  ratio  of  0.05  and  Reynolds  number,  based 
on  the  nose  radius,  of  Re  = 430. 

03 

Initially  the  shock  was  assumed  parallel  to  the  body  surface 
and  at  a distance  of  0.1072.  This  value  is  the  converged  stagna- 
tion point  shock  stand-off  distance  obtained  using  Davis'  method. 


28 


A ED  C TR-79-25 


The  initial  shock  shape  is  in  general  not  critical  to  the  overall 
solution  and  other  smooth  initial  shock  shapes  could  be  used 
as  well.  With  the  initial  shock  shape  defined,  local  iterations 
were  necessary  to  solve  the  star  sweep  equations  for  the  new 
shock  shape.  Figure  (2)  shows  the  normal  velocity  component,  v, 
at  n = 0.9  9 at  several  s-locations  for  each  iteration,  v/here  it 
is  observed  that  the  solution  converges  at  all  locations  in  less 
than  six  iterations.  The  solution  did  not  only  converge  faster 
than  the  solution  in  reference  [5] , but  also  was  more  stable 
since  the  oscillatory  behavior  observed  in  the  same  reference 
was  completely  eliminated. 

The  choice  of  the  time  step  size,  in  the  alternating 
direction  method,  is  analytically  a hard  task.  In  such  cases 
numerical  studies  can  be  conducted  to  find  the  basis  for  the 
choice  of  the  time  step  size.  Figure  (3)  presents  the  results 
of  a numerical  study  conducted  to  show  the  convergence  of  the 
shock  radius  with  the  global  time  iteration  cycles  for  different 
time  step  sizes  ranging  from  it  * 10  to  At  = 60.  If  the  solution 
is  assumed  converged  when  the  change  of  the  shock  radius  does 
not  exceed  0.1%  of  its  value,  the  figure  indicates  that  this 
convergence  limit  is  achieved  faster  with  the  largest  time  step 
size.  At  = 60.  On  the  other  hand,  if  the  convergence  limit  is 
set  to  be  0.07%,  the  results  show  that  this  limit  is  achieved 
faster  with  the  smallest  time  step  size,  At  = 10.  Therefore, 
the  choice  of  the  time  step  size  will  generally  depend  on  the 
degree  of  accuracy  desired. 
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In  addition  to  the  factor  mentioned  above,  the  consistency 
limitation  of  the  alternating  direction  implicit  method  provides 
another  factor  in  the  time  step  size  choice.  The  consistency 
limitation,  as  described  in  reference  [5],  shows  that  the  two 
half  time  sweeps  converge  to  slightly  different  values. 

Figures  (4)  and  (5)  represent  the  converged  shock  radius  in 
the  two  sweeps  for  two  different  time  step  size  values,  namely 
At  = 10  and  At  = 60,  and  with  a convergence  limit  of  0.1%.  If 
we  compare  the  two  figures,  we  find  that  the  difference  between 
the  two  sweeps  is  more  noticeable  when  large  At  values  are  used. 
Therefore,  for  this  difference  to  be  minimum,  small  time  step 
sizes  are  always  preferred. 

On  the  other  hand,  a lower  limit  was  observed  on  the  time 
step  size  beyond  which  the  solution  diverged  during  the  global 
iterations.  This  limit  depends  on  the  initial  guess  of  the  shock 
shape  and  on  how  far  the  calculations  are  to  proceed  downstream 
on  the  body  surface.  As  a conclusion,  it  is  clear  that  the 
proper  time  step  size  for  the  present  method  is  the  smallest 
value  which  yields  to  a converged  solution. 

Comparisons  of  the  present  results  with  Davis  and  Werle, 
et  al . , references  [4]  and  [5],  are  shown  in  figures  (6)  to 
(8).  The  calculated  wall  pressure  distribution  in  figure  (6) 
shows  that  the  present  results  agree  well  with  the  other  cal- 
culations. Figures  (7)  and  (8)  also  indicate  that  good 
results  for  the  surface  skin  friction  coefficient  and  the 
Stanton  number  are  obtained.  In  general,  the  agreement  between 
the  results  of  the  different  methods  are  remarkably  good  and 
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they  essentially  reproduce  one  another. 

To  complete  the  comparison,  another  test  condition,  with 
a high  Reynolds  number,  was  selected.  Figure  (9)  shows  the 
calculated  results,  together  with  the  experimental  data  of 
Little  [12],  for  a flow  over  10°  half  angle  hyperboloid  with 
a Reynolds  number  of  63,800,  free  stream  Mach  number  of  10.12, 
wall  to  stagnation  temperature  ratio  of  0.9,  and  a free  stream 
temperature  of  90 °R.  The  comparisons  are  remarkably  good 
ever  the  entire  body  length,  specially  downstream  at  a distance 
three  times  the  nose  radius. 

Based  on  all  the  previous  comparisons,  the  present  method 
is  comparable  in  accuracy  with  the  other  methods  and  is  capable 
of  predicting  excellent  results  without  the  oscillatory  behavior 
observed  in  reference  [5]  and  without  the  stability  limitation 
on  the  number  of  the  global  iterations  observed  in  reference  [10] . 
Few  other  cases  were  also  investigated  to  show  the  applicability 
of  the  method  to  regions  far  away  from  the  stagnation  point 
and  the  capability  of  handling  slender  body  configurations. 

Figures  (10)  , (11)  and  (12)  show  the  pressure,  skin  friction 
and  Stanton  number  distributions  over  10°  and  22.5°  half  angle 
hyperboloids.  The  results  of  these  test  cases  were  obtained 
using  a time  step  size  of  fit  = 600  and  400,  and  a spatial  step  size 
of  As  = 0.1  and  Ag  = 0.02.  The  numerical  calculations  were 
performed  on  Amdahl  470  V/6  computer,  and  the  computation  time 
was  25  minutes.  With  these  results,  it  is  clear  that  the  method 
was  not  only  very  successful  in  extending  the  region  of  cal- 
culations far  downstream  without  any  numerical  difficulties,  but 
also  successful  in  handling  slender  body  configurations. 
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VI.  CONCLUSIONS 

A preliminary  investigation  has  been  made  of  the  use  of 
a coupling  scheme  to  solve  the  fully  viscous  shock  layer 
equations.  The  time  relaxation  technique,  of  reference  [11], 
was  also  incorporated  into  the  scheme  to  allow  for  the  shock 
shape  changes  during  the  iteration  process.  Through  application 
of  the  method  to  hypersonic  flow  past  hyperboloids,  it  was 
found  that  the  method  is  stable  and  produces  results  in  good 
comparison  with  other  methods. 

The  main  success  of  the  scheme  lies  in  its  ability  of 
extending  the  flow  field  calculations  far  aft  of  the  nose 
region  without  the  difficulties  observed  with  the  cascading 
scheme  of  references  [4]  and  [5],  The  method  is  also 
capable  of  handling  slender  body  configurations  where  a 
large  portion  of  the  flow  field  in  the  shock  layer  will  be 
predominantly  inviscid. 

There  is  no  difficulty  in  the  choice  of  the  shock  initial 
shape  and  shapes  parallel  to  the  body  surface  with  constant 
thickness  are  very  adequate.  However,  the  time  step  size  is 
of  fundamental  importance  to  the  solution  convergence  and  it 
depends  mainly  on  how  far  downstream  the  flow  field  is  to  be 
investigated. 

Difficulties  were  encountered  when  the  method  was  applied 
to  sphere  configurations  and  hyperboloid  configurations  with 
high  Reynolds  number.  In  these  cases,  variable  grid  size  across 
the  shock  layer  must  be  used  and  further  studies  are  needed  for 
the  proper  choice. 
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APPENDIX  A 

v 

SHOCK  DERIVATIVES 


The  shock  derivatives  dus^/ds,  dT^/ds,  dpg^/ds,  and 
dv^/ds  have  been  presented  in  detail  in  Appendix  A of 
reference  [5]  and  their  derivation  will  not  be  repeated.  For 
present  purposes,  the  results  are  only  summarized  and  are 
given  by: 


dUsh  _ da  d^_ 

ds  1 ds  2 ds 


(Al) 


^Psh  „ da 
ds  3 ds 


(A2) 


dT  . , 

sh  „ da_ 

ds  *4  ds 


(A3) 


and 

dVsh  _ da  dp 

ds  " R5  ds  6 ds 


( A4 ) 


where 


y — 1 <=]n  v — 1 1 

K.  = (d— — - 1)  'sin(2a  - <p)  + sina  sin{a  - *)  - — K-. 

1 2 psh  7 1 

• y-i  ^sh 

K_  = cosa  sm  {a  - p)  - sina  cos  (a  - <f>) 

2 Y psh 


K,  = 


3 (Y+l) 


sin  2a 
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2y  4 cosa 

K.  = ! — ~ sin  2a  + s— - ^ 


(Y+l) 


(Y+l)  HZ  sin  a 


T i 

K-  = (1  - r^)  cos  (2a  - 4>)  - sina  cos(a-4>)  ^T1  Ki 

3 Y Psh  ' 


and 


K,  = -cosa  cos(a-f)  - sina  sin(a  - <p) 


T 

Y-l  sh 


Y P 


sh 


where 


K1  = 


2yV 

Tr^Ty+TT 


sin  2a  + 


_4y cosa 


3 4 .2 

4y  M sin  2a  sm  a 


(y  + 1)M^R^  sin^a  R2  (Y+l) 


2yM  sin  2a 


Y (Y+l)  _ 2y  (Y-h  1 , 4y  sin  2a 

Y-l  Y + l v ~2  _2 


(Y+l)  sin  a 


and 


= [2y  M2  sin2a  - (Y~l)] 
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APPENDIX  B 

SHOCK  SLOPE  DERIVATIVE 


The  shock  slope  derivative  is  derived  in  this  appendix, 
for  use,  together  with  shock  derivatives  of  Appendix  A,  in  the 
viscous  shock  layer  solution. 

In  the  spatial  coordinate  system  the  shock  angle,  a,  is 
written  as,  (Figure  1) 

a = tan  d (4^ — ) (Bl) 

dxsh 


where  R = yB  + n^cos#  and  xgh  = ^ - nghsin? 


( B 2 ) 


Hence  the  derivative  da/ds  is  evaluated  as. 


da 

ds 


,2_  dx  . 
d R sh 


Note  that 


dxsh 

5i cos+(l+Kngh)  -nsh  sin<t» 


and 


dnsh 

— = (1+Kngh)  tan(a-<f>) 


combining  ( B 4 ) and  (B5)  yields 


dx 


sh 


ds 


<1+,<nsh) 


cosa 
cos  (a-$) 


(B3) 


(B4 ) 


(B5) 


(B6) 
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2 2 

Hence  the  derivative  d xg^/ds  is  evaluated  as, 


,2 

d x 


ds 


Sh  = (K,n  ■+Kn*  ) + {1+Kn  ) [.  gintt  p- 

2 sh  sh  cos  fa-$ ) sh  cos  (a-tf ) ds 


cosa 


cos  ( a — 4>  > 


sin  (a-<J>)  (||  + k)  ] 


(B7) 


substituting  (B6>  in  (B3)  and  noting  that  dR/dx^  = tana 
yields  after  certain  manipulations, 


da  _ ,,  * cos^a 

ds  Knsiv  cos(a-<{>) 


d2R 


(B8) 


dx 


sh 


It  is  to  be  also  noted  that. 


tana  = 


dR 


dR/ds 


^sh  dl£sh/ds 


<B9) 


Hence  the  second  derivative,  d^R/dx2^,  can  be  shown  to  be 


d2R  _ d2R/ds2 


dx 


sh 


(dxgh/ds) 


d2x  . /ds2  dR/ds 
sn 

(dxsh/ds)3 


(BIO) 


2 2 

Substituting  for  dxgh/ds  from  (B6) , d xeh/ds  from  (B7) , 


‘sh7 


and  dR/ds  from 


dR  - n 4-k^n  i sing 
ds  s cosTa-TT 


and  then  evaluating  (B8)  yields  after  proper  manipulations. 


2 2 

da  _ d R r cos  (a-<5)  ■■  _ dR  ,icsin  ( 2ot —2 ^ ) + 

ds  ^g2  (l+KnoV,)  cos(j)J  ds  lcas<j>  { l+icn,,^) 


.2 

> 

lsh 


nsh<  'cos 

j-] 


sh'  cos<J>  (l+krng^) 

(Bll) 
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APPENDIX  C 

ADI  FORMULATION  OF  S -MOMENTUM  EQUATION 


The  s -momentum  equation  is  written  in  the  following  nor- 
malized form 


ifu 

3n2 


+ a. 


3U 
3 n 


+ a. 


u + a. 


+ 


3u  _ 
Q rr  — 0 
4 3 £ 


(Cl) 


where  a^,  a^,  and  are  defined  by  equations  (5a-d) . Note 

• i 

that  ush  and  pgh  appear  in  the  coefficients  and  a ^ . Upon 
using  equations  (Al) , (A2)  and  (Bll)  into  equation  (Cl)  and 
after  certain  manipulations  the  following  equation  is  achieved, 
(details  are  given  in  reference  11) , 


92u 


a , B O 

61  ^ “2  ds2 


■ dR 

'3  ds  34 


8 1H  - n 
35  3£  ° 


(C2) 


where 


Y2u  + 


Y3u  + 


Y4u  + 


Y 

Y 

Y 


5 

6 
7 


and 


Y 2 > Y 3 r y4.  y5>  Y6  and  Y?  are  given  by 


V = -A.  K COS  (a~^}  - 

t2  *1  (l+icn  . ) cos  4 


sh 


• /n  n . K'cos  (a-9) 

Y = +A  K r*sin(2g-2»)  + _sh 

y3  A 1 lcos<t  (l+Kn  , ) ^ 


sh  cos  if)  ( l+<ns^) 


(C3) 
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Y4  = +Aic  + B + C + D 


Tc  = 


A P K_  cos  ( Ct— (J> ) 


1 psh  3 (1+KlWcos 


sh' 


_ , P v r <sin ( 2«-2$)  s 

1 ^ ^ ^ A / 1 i . \ 


n k 1 cos 


1 psh  3 co^(l+.cnsh>  cos^1+<nsh)2 


y7  = A (p  - — ^ n p ) 
7 1 5 nsh  n 


where 


A = 


p . n , <n  . 
sh  sh  sh  pu 

2 l+<n  , ri  - 

e Psh  sh  u 


B = - 


p shvshnsh 


2 

c y 


sh 


II  , 

sh  p v 

l+Kn  *,1  - 
sh  y 


>cn 


C = - 


sh  yn 


^^sh71  y 


icn 


D - - (- 


sh 


cos  4 n 


sh 


<n 


sh 


1+Knshn)  r+nshcos<j>n  (l+<nshn) 


. _ Pshnsh  nsh 

1 ~ 2 


(1+kh  , n)  - 
sh  sh  w Ush 


It  can  be  shown  by  proper  substitution  that 


n ,<* 

r2  - - [2K  tan(a-*>  + 

2 sh 


(C4) 


(C5) 


(C6) 


52 


A E DC-TR  -7  9-25 


APPENDIX  D 

GOVERNING  EQUATIONS  LINEARIZATION 


s-Momervtum  Equation : 


2- 

3 u , , 3u.  , - , , / 3 u> 

7"*2  1 W + a2u  + a3  + a4  {IT) 

3 n 


= 0 


(Dl) 


where  a and  are  given  by  equations  (5a-d) , and  they 

are  rewritten  as, 


= a11  pu  + a12  pv  + a13 


a „ = 


°21  °u  + “22  PV  + “23 


n 


a = 


a 31  (PC  " n 


sh 


n P„  + 


sh  - 


sh 


n p 


P) 


sh 


a4  = a41  PU 


(D2) 


and 


°shushnsh  nsh 


11 


e ii 


sh 


1+icnshr| 


pshvshnsh  1 


12 


£ JJ 


sh 


_ T1 


<n 


sh 


cos<t>  n 


sh 


13 


- l+<nshn  r+nsh  ncos< 
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pshnsh 
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nsh  ush 
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sh 
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p shvshnsh  nsh  1 
2 l+<n  . - 


a22 

2 

£ vsh 

1+,:nsh 

M 

nsh 

V 

nsh 

cos<|>  n . 

I Sh 

°23 

< 1+Knshn 

1+Knshn 

r+nshr'COS^' 

P3hnsh  11 

sh  1 
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a 31 

2 1+icn  , n 

e ^Sh  sh  v 

ush 

p shushnsh 

nsh 

a41 

2 

e ^sh 

1+Knshn 

P 

Using  the 

linearization 

formulas  given 

by  equations 

equation 

(Dl)  could  be 

written 

as. 

■) 


(22a-c) 


(D3) 


U + a 


o J_u> 

4 sr 


+ ( 


,3uu 
1 3n 

L -o 

■ + u <*2  + 

3u 

35 

a4) 

/ o 
(al 

3u°  . o-o, 
Sn  + °2  u + 

O 

“4 

3U°. 

35 

= 0 


(D4) 


where  the  quantities  with  superscript  ° are  evaluated  from  either 
the  previous  iteration,  or  the  previous  5 step.  Upon  using 
equations  (D2)  and  (D3) , substituting  for  the  density  from  the 
equation  of  state  p = pt,  and  using  the  linearization  formulas 
(22a-c) , equation  (D4)  can  be  written  in  the  following  form, 


7 T + U1  + a3  V + (bl  "n  + b3  V 
3 n 

+ (c^  u + c2  v + p + t)  + d = 0 (D5) 
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where  the  coefficients  a^„  a^,  . etc.,  are  given  by  the 
following  f orms r 
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al  = " Y1Y2  v 


sh  1 -o  -o 
— P u 
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sh  sh  sh  y 


u , n ' 

, , sh  sh  n -o  -o 

b,  = y,  lv,  - — - — — p u 

Sh  sh  y 


-o  -o 

3 V 


] + IT1  + <V2  + Y3 

y u 


, _ , sh “sh  n 

3 Tir2  pshushnsh  nsh  y 


r -o  Ush  u°  t v°  Ush  nsh  r,  3u.0 

C1  = - Y2  [y1  p <2  — - - — + < — - — - — ~ — 

sh  y y sh  sh  y 

u . , _-o  y 

sh  l 3u  . , , n . s i 

+ VT  ~ W]  + K{~  + Ky2  + Y 3 
sh  y y 


-o,l  ou° 

c2  ( _ gn 

y 


-o 

+ <y2  r-) 

y 


c,  = - 


-o  u 

Y1  |o  Iy2<F 

t £ 


n u 

sh 

n 

5U° 

ush 

n v 
sh 

y 

•s  „ 

a r 

vsh 

1 -o 

— u 


sh  1 


sh  y 


-o  . „-o 

V . -O  1 oU  -o, 

< — u - — -s v i 

- -or, 

y y 


4 ^sh  1 

c3  -o  “ Ylr2  p u . v . - 

o sh  sn  sn  y 


C I 


AE  DC-TR  -79-25 


and 


, . “sh  “sh  35°  , ush  -o  , “sh  a!°  S°  u° 
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sh  “sh 
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(D6) 


where 


Y1  = 


pshVshnsh 


£ P 


sh 


n 


sh 


1+Knshn 


Y i = 


n . cost 
sh 

r+nsh  ^,CCIS', 


(D7) 


Energy  Equation: 

The  energy  equation  is  given  by  equation  (6)  and  is  written 
in  the  form, 


5 2 1 
9n2 


3 1 
3 n 


+ a_  t +•  a,  + a. 


at 

95 


= 0 


(D8) 


where  the  coefficients  ct^,  and  oc^  are  given  by  equations 

(6a-d) , and  they  are  rewritten  as. 
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1 

2 


an 

a21 


pu  + 

pu 


12 


pv  + a 


13 
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and 


a 3 

n v 

sh 

<p< ' 

- Psh 

n p + - — 

n psh 

p)  u + a32  v pn 

+ “33  (S„ 

*nsh 

u>  2 

1+Knshn 

II 

«41  P u 

1 

p , u . n , 0 
sh  sh  sh 

nsh 

n 

°11 

2 

E psh 

1+Knsh 

V 

“12 

pshvshnsha 

_1_ 

2 

£ ysh 

p 

(D9) 


13 


_ n 


xn 


. cos^n  , 1 

sh  sh 


- l+<nghn  r+nshncoS* 


pshushq 


n 


sh 


sh 


21  2 _ l+<n  . n 

' "shTsh  sh  " 


PshUshnshq  nsh  1 
“31  2 _ l+Kn  . n 

e vshTSh  sh  11 


P iV  un  uff  1 

rsh  sh  sh  1 

‘32  ,2„  t 7 

e yshTsh  u 


2 

ushq 


33  T 


sh 


p . u , n .a 
sh  sh  sh 


n 


sh 


‘41 


£ U 


sh 


l+<n  7 n - 
sh  v 


(DIO) 
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The  linearization  of  equation  (D8) , results  in  the  following 
equation 


ift 

3n2 


, , o 3t  , o - . o 3t. 
+ (“i  t—  + a-  t + a,  — ) 


1 3n 


4 H' 


+ ‘if  ‘l  + _t°  »2  + °3  + if  °4> 


, O 3 1°  . o-o,  o 3t°.  _ n 

(c,i  — + “2 1 4 "4  r'  0 


(Dll) 


Upon  using  equations  (D9)  and  (DIO),  substituting  for  the  density 
from  the  equation  of  state  p » pt,  and  using  the  linearization 
formulas  (22a-c) , equation  (Dll)  can  be  written  in  the  following 
form 


32t 


— | + (a3  P5  + a4  t^)  + (bj^  un  + b3  Pn  + b4  tn) 
3 n 

+ (c,  u + C-  V + c,  P + c.  t)  + d = 0 


(D12) 


where  the  coefficients  a3,  a4,...  etc.,  are  given  below 


a-.  = y,y 


Ush  ^sh  _o_  -o 


1'2  v*.  P 


sh  sh  sh  v 


a,  = - 


bl  " 2 


b3  = Y- 


Y1^2 

ush 

vsh 

a -o  -o 

— P u 

V 

2 

ush° 

(u° 

T1 

-o. 

- <Y2  u ) 

Tsh 

a 

^sh 

( Ush 

u 

pshT 

sh  vsh 

n 


2 n 


sh  n u°  + v°) 
sh 
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u , n , 

, , sh  sh  o -o  -o 

b4  = Y1  v7”  y2  — n D u 

sh  sh  y 


a -o-o. 
— P v ) 


+ r21  + <y,  + y. 


C1  = y2  °{y1[p  * 


„ u . n . 
-o,  sh  sh 


v . n , 
sh  sh 


3 1°  Ush  Tsh  -o 


3 n 


v . T , 
sh  sh 


a , ,r3  . 

sh  3 t . 1 

v .35  - 

sh  y 


Ush  Psh  1 , — o nsh  -o 

~ J— T“  T ,PE  ' 5-  n pn 


sh  sh  sh  y 


sh 


' psh  -o.  . - Ush  .,-o  -o.  , 

+ - — p ) ] - 2k  - — (u  - k y-  u ) } 

psh  Tsh  n i 


c,  = r 


(_  if.  3tf.  + psh 
1 ; t°  5"  »shTsh 


?n> 


a r ~or (ush  nsh  3t°  Ush  ^sh  zo 

c,  = y,  — iy,  u [(- — - — n -r- = — t 

3 1 y 2 VSh  nsh  ^ 


vsh  Tsh 


Sh3t°,  1 

vsh 


u 


sh  Fsh 


v°  3t° 


j.  an  an  i vi-  | 

1°  vsh  pshTsh  ' t°  9x1 


C4 


„ „ u . n , -o  u . ,to 

-o-o.  sh  sh  at  sh  3t  , 

■ h-  ‘*2  ’ “ <7^  » — - 7^  5T1 


— o — o at  i l 
- p V -TT-i 


3n 


y t° 


and 


u . n . -o  u , T , 
a , r-o  -o  , sh  sh  at  sh  sh  -o 

d = - y,  — (y,  [o  u (- — - — n 7- — s? — t 

1 y 2 VSh  nsh  °n  VSh  TSh 


^sh  + %h  Psh  -o  (-o  ^sh  n po  + 

vsh  ^ vsh  shTsh  ? nsh  n psh 


sh  -o 


P )]  " 
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- ?°  If  + v°  p°}  - ^r1  ^ “0) 


sh  sh 


sh 


where  y^ , y2  and  Y-j  are  given  by  equations  (D7)  . 


n-Momentum  Equation: 

The  normal  momentum  equation  is  given  by  equation  (8) 
rewritten  in  the  following  form 

1 I 

e (-  ^ n v + Is**  v)  + e v 

1+Knsh  5 nsh  11  vsh  ush  nsh  71 


* Ugh  P u2  + - Psh- 


1+‘nshn  vsh 


pshushvshnsh  n 


where  with  full  shock  layer  6=1,  and  with  the  thin  shock 
approximation  6=0.  Using  the  equation  of  state,  p = pt, 
equation  (D14)  can  be  written  in  the  following  form 


(1+Kng^n)  pshushvshnsh  ^5  ngh  n Vn  + vg^ 


-EE. 


2 _ — — - 2 < --2 

+ P . v . 0 p v v — p , u . n . •=— p u 

sh  sh  v n sh  sh  sh  l+<ng^ri  ^ 


+ p , t p =0 

^sh 


The  linearization  of  this  equation  gives. 


c u + (a-  v + b v + c_  v)  + (b  p + c p) 
1 2 £ 2 7i  i jnJ 


+ (c4  t)  + d = 0 


(D13) 

and  is 

(D14) 

layer 

(D15) 

(D16) 
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The  coefficients  c,  , a,,,  etc.,  are  given  by  the  following 

relations , 

f 

n , 

c,  = y,p„uu  . [-2<u  , p°  u°  + s v , (v?  - — — n v° 

1 2 sh  sh  sh  r sh  £ n , n 

, sh 


+ v°)  p°] 

vsh 


a2  = 7 2 pshVshUsh0  P° 


-o 

n u. 

, sh 

-o 

b2 

pshVsh  6 

P 

{-Y-  n u 

2 nsh  sh 

u + v , 
sh 

-o 

i 

v , 

, sh  -o 

-o 

C2 

pshVsh  5 

p 

Yn  — U , U 

1 2 v . sh 

sh 

+ V . V 

sh  n 

-o, 


b,  = 


psh  *° 


^ n , v , 

t r />  sn  — O SI1  wOi  — o 

c = P .v  {y  [0  u . (v  - - — n v + - — v ) u 
J sh  sh  2 sh  i n . n v. 


sh 


sh 


- K U°2  ] + 9 V v°  v°  } 

” sh 


sh 


psh  pn 


and 


n 


d . = 


sh 


psh  Pn  + 2 pshVsh  ^l+|'ns[1n  “sh 


-o  -o  r sh  -o 
P u - — u 


sh 


- 0 (v°  - — ^ n v°  + v°)  ] - 0 v p°  v°  v°}  (D17) 

5 nsh  11  vsh  sh  n 


where  ■'■s  gi-ven  by  equation  (D7)  . 
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Continuity  Equation; 

The  continuity  equation  is  given  by  equation  (7)  and  after 
certain  manipulations  it  can  be  written  in  the  following  form 

P(a1  ^ + “2  ^ + a3  + a4  % + “5  V + a3  5 pc 


+ (a4  v + a5  u)  = 0 


<D18) 


Upon  using  the  equation  of  state  (9) , equation  <D18)  can  now  be 
written  as 

pt  {a^U  + ct2  V + a3  U^  + a4  + a5  u^) 


+ a3  u(t  p^  - P t?)  + (a4  v + a5  u) (t  - p tn>  = 0 


(D19) 


where  the  coefficients  a^,  etc.,  are  given  below, 

tnshE-"sh  ujh  P,h  (r+n  h„co*») 

U_,_  [ — —* - — + - — + T^rrz — ^ 


a,  - 


sh 


1 1+<nshr'  sh  sh 


sh 


ush  psh  (r+n^ncos*) 


n , ncos<ti 
sn 

(r+nstincosi(>) 


and 


= 


a , - 


a „ = 


n , cos<j>  icn  , 

r sh  __sh_, 

pshvsh  (r+nghnco"s4>y  1+|i:nsh 


“sh 

pshush  1+icn^n 

p .v  , 
sh  sh 

n , n , 
sh  sh 

psh  ush  nsh  n l+icns^n 


(D20) 
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where  the  shock  derivatives  n , and  n , will  result  from  the 

sh  sh 

first  and  last  terms  of  equation  (8) . In  the  present  calcula- 
tions, 'n  . was  evaluated  as  the  shock  thickness  derivative  at 

Sh  y. 

^ , 1 n 

the  time  t*,  while  n . was  evaluated  at  the  time  t . "In  the 

sh 

limit,  when  the  solution  convergence  is  achieved,  ng^  value 
T 5 

approaches  n ^ value. 

The  linearization  of  equation  (D19) , using  equations 
(22a-c),  gives  the  following  equation. 


(a,u  + a^p.  + a,t  ) + (b,u  + b_v  + b,p  + b,t  ) 
1 £ 3 1 4 5 In  2 n 2 t\  4 n 


+ (c^u  + c^v  + c^p  + c4t)  + d = 0 


(D21) 


where  the  coefficients  of  the  linearized  equation  are  given 
below, 

al  = y2  p sh  ush  ^ b 
a 3 = Y 2 psh  ush  “°  S° 


-o  -o 

a . = -y_  p . u . u p 
4 2 sh  sh  r 


n 

■ Sh  ,,  rO  -o 

b,  -Y->  P u . ti  p t 

1 2 nsh  sh  sh  r 


, -o  -o 

b2  = psh  Vsh  P 1 


n , 

, , -o  sh  -o.  -o 

b3  = (psh  vsh  v * psh  ush  ^ n y2  u ) t 


h , = 


- b3  p°/t° 
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C1  = psh  ush  y2  {[— 


* v,  <- 


sh  ~nsh 

k Ush  + 

psh 

nsh 

ush 

psh 

1 

COS  <j> ) 

- i 

n)  ] 

cos  4 

nsh 

- „ (t°  5°  - i°  t°)  ♦ Ct°  p°  - p°  t?) ) 


c2  ■ »sh  vsh  < <‘t0  5;  - 5°  E°>  + '^3  + «T2)  P t ! 


»sh  vsh  & + E°  - E°] 


+Y,  ( 


"sh  y2  tE° 

<■**» 

ci 

0 

1 

n 

sh 

nsh 

u°) 

n 

nshr-nsh 

S°  [ „5 

nsh 

f 

+ Ush 

+ — 

ush 

p sh 

1 

{r+nshncos4) 

nsh 

nil  - 

t 

nsh 

nshcos4 

nsh 

Uf  ] 

nsh 

-2  »5h  v3h  [5°  t°  5°  + t«72+73)  5°  5°  t°  * v°(t°  5°  -5°  t°)l 


-2  f.h  ush  y2  (P°  - n + “°(E°  5°-  EC> 


-c  -o  -o  "-t-"-  /5ht!sht,  (r^cost)  _ngj, 
P ' "sh  "sh  psh  3 "sh"05*  "sh 


- -£±L  n U°  (t°  P°  - 5°  t°)  } 

nsh 
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APPENDIX  E 

FINITE  DIFFERENCE  EQUATIONS 


This  appendix  gives  the  coefficients  of  the  finite  difference 
equations. 

s-Momentum: 

Using  the  difference  quotients  (27-29),  equation  (23)  can  be 

/ 

written  in  the  following  difference  form. 


{ a , , u . + 

11  m,n-l 


a.  ,p  ,)  + (b..u  + b.  _v  + b.  ,p  „ 

13  m,n-l  11  m,n  12  m,n  13  m,n 


b.  .t  ] 
14  m ,n 


+ (c.-u  .1  + c . ,p  „ , . ) = d, 

ll  m,n+l  13^m,n+l  1 


After  dropping  the  subscript  m,  this  equation  becomes: 

,allun-l  + *13P„-l,+(bllun  + b12vn  + b13pn  + b14V 


+ (cuVl  + °13pn+l)  ' dl 


(El) 


where 


a. , = (2  - b.  in  ) /An  , (An  + An  .) 
11  1 n n-1  n n-1 


i.  , = - b-  An_/An_  , (An_  + An_  ,) 
13  3 n n— l n n-J. 


b . = c.  + [-2  + b,  (An  - An_  .)]/An  An  , + a./(A£  CRNI) 
11  1 In  n-1  n n-1  I 


b12  ' c2 

b13  ’ a3/U5  CHNI>  + c3  * b3 


b.  . = c , 
14  4 
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c, , = (2  + b.  An  i ) / An  (An  +An  .) 
ii  i n-i  n n n-i 


C13  = (b3  Ann-l)/inn  (Ann  + A71n-1} 


1,  = R + (a.u  . + a,  pm  J/Ue  CRNI) 

1 1 m- 1 , n 3 m- 1 , n 


(E2) 


The  coefficient  CRNI  is  equal  to  0-5  for  Crank-Nicolson  scheme, 
and  1 for  pure  implicit  schemes.  The  coefficient  R is  given  by, 

R = - (d/CRNI)  - [(l-CRNI)/CRNIl[(unn)iti_^n  + 

+ C1  um-l , n + c2  vm- 1 , n + c3  pm-l,n 
+ c4  ^-1,  n 


Energy: 

Using  the  difference  quotients  (27-29) , and  dropping  the 
subscript  m,  equation  (24)  can  be  written  in  the  following 
finite  difference  form. 


(a21un-l  + a23pn-l  + a24Vl>  + (b21u, 


b22vn  + 


b23pn  + 


b24tn> 


(c21un+l  * c23pn+l  + c24tntl) 


(E3) 


where 


21 

' bl 

An_/An_  , 
n n- 1 

Unn  + 

Ann-1} 

23  “ 

- b3 

An/An  . 

n n-i 

(Ann  + 

Ann-1] 
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a24  = (2  ' b4  Ann)/A"n-1  Unn  + Ar,n-1] 


b21  * C1  + bl  Unn  " Ann-1)/Ann  Ann-1 


b22  C2 


b23  = c3  + a3/A?  CRNI  + b3  (Ann  * 'Vl*  AVnn-l 


b24  = C4  + a4/A?  CRNI  + b4  (Ann  " Ann-1) /AnnAT1n-l 

- 2/in  An  , 
n n-1 


1,  = R + (a^p  , + at  _)/A? 

2.  j m-i,n  4 m-1 , n 


CRN! 


and 


R = - 


d/CRNI  - [ ( 1-CRNI ) /CRNI ] [ ( t ) , „ + b. (u  > - r 

nn  m-l  ,n  l ti  m-l,r 


+ b_(p  ) , + b . (t  ) . + c.u  , 

3 n m-l,n  4 n m-l,n  1 m-l,n 


+ c~v  . + c,p  . + c.t  , ] 

2 m-l,n  3rm-l,n  4 m-l,n 


n-Momentum: 

Using  the  difference  quotients  (30,  31),  evaluated  at 
(M,n+^) , and  dropping  the  subscript  m,  equation  (25)  can  be 
written  in  the  following  form, 


(b,,u  t b,_v  + b,,p  + b,.t  ) 
31  n 32  n 33rn  34  n 


+ (c,tu  ,,  + c--v  + c.-p  ..  + c_.t  . , ) = d_ 
31  n+1  32  n+1  33  n+1  34  n+1  3 


where 

b^  = CRNI  A£ 


(E4) 


(E5 ) 
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b32  = a2  + CRNI  AC  c2  " 2 CRNI  b2 

b33  = CRNI  4?  Cj  - 2 CRNI  t>3  AC/An 

b34  = CRNI  AC  C4 

C31  - CRNI  AC 

c32  = a2  + CRNI  c2  + 2 CRNI  b2  A?//An 

c33  = CRNI  AC  c3  + 2 CRNI  b3  AC/An 
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where  the  coefficients  a^,  a42'  ebc  • • are  ,?iven  asr 

a . , = a,  + CRNI  AC  c,  - 2 CRNI  bn  AC/An 

41  l j.  J- 

a . _ = CRNI  4?  c,  - 2 CRNI  b,  AC/An 

42  2 2 

a,,  = a,  + CRNI  45  c,  - 2 CRNI  bQ  AC/An 

43  3 3 J 

a.  , = a.  + CRNI  AC  c,  - 2 CRNI  b.  AC/An 

4 4 4 4 4 

b41  * aL  + CRNI  AC  cL  + 2 ,CRNI  ^ A C/An 

b42  = CRNI  AC  c2  + 2 CRNI  b2  AC/An 
b43  = a3  + CRNI  AC  c3  + 2 CRNI  b3  AC/An 

b44  = a4  + CRNI  c4  + 2 CRNI  b4  A^/An 

d4  = - ! H a * [il-ola-cMm1i]|u,.lrnHtVi,n) 

- o2(1-CRNI)4£  (vn.lin+1  + vm_1/n> 

+ [«3  - =3I1-CFHI)45]  <pm.ljtl+1  + Pm.1/n> 

+ [»4  - c4(l-CRBI)15l  (tm.lin+1  + Sn-l.n’ 

- 2 (1-CRMI)  fi  n-iVl.w!  * um-l,n> 

+ b2(vm-lfn+l  “ vm-l,n}  + b3(pm-l,n+l  " pm-l,nJ 
+ b4(tm-l,n+l  " ‘m-l.n11 


(E8 ) 
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APPENDIX  F 

SOLUTION  OF  THE  DIFFERENCE  EQUATIONS 

The  difference  equations  (32)  to  (35)  are  arranged  into  a 
system  suitable  to  inversion  with  the  recursion  formulas  (38a-d) . 
The  boundary  conditions  (36a-c)  and  the  n-momentum  equation  (34) 
evaluated  at  n=l  are  written  together  to  obtain. 


u^  = 0 


v,  = 0 


t = t 
1 w 


(b31ul  + b32vl  + b33pl  + b34bl> 


+ (c3J.U2  + C32V2  + =33P2  + c34b2> 


= d_ 


(FI) 


While  the  boundary  conditions  (37a),  (37c)  and  (37d) , and  the 
continuity  equation  (35)  evaluated  at  n=N  are  written  together  to 
obtain 


"M  " 1 


<a41UN-l  + a42VN-l  + a43PN-l  + a44tN-ll 

+<b41  “n  + b42  VN  + b3PH  + b44tN»  = d4 


At  the  shock  n=N,  we  still  have  an  additional  boundary  condition, 

v„  = 1.  The  value  of  the  shock  thickness  i's  iterated  on  to 
N 

/ 
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satisfy  this  boundary  condition.  The  finite  difference  equations 
(32)  to  (35)  evaluated  at  n = 2,  3,...  and  (N-l) , are  given  by 

(allun-l  + a13Pn-l}  + (bllun  + b12vn  + b13Pn  + b14V 


+ (cllUn+l  + C13Pn+l)  " dl 


(F3) 


(a21un-l  + a23Pn-l  + a24tn-l)  + (b21un  + b22vn  + b23Pn  + b24tn) 


+ {c21Un+l  + C23Pn+l  + C24tn+1)  + 


(F4) 


ib31Un  + b32Vn  + b33Pn  + b34V 


+ (c31Un+l  + C32Vn+l  + C33Pn+l  + C34tn+1)  = d3  (F5) 


(a41Un-l  + a42vn-l  + a43pn-l  + a44tn-l) 


+ (b41Un  + b42Vn  + b43pn  + b44tn)  = ^ 


(F6) 


Equations  (FI)  to  (F6)  in  this  arranged  form,  represent  a system 
of  block  tridiagonal  equations.  Their  solution  is  obtained  with 
the  following  recursion  relations, 

= D u+E  v+F  p+G  t+H 
un+l  un+l  n un+l  n Un+1  n un+l  n Un+1 


V 


„ = D u+E  v+F  p+G  t+H 
n+1  v„^  n vn+1  n vwi1  n n v 


'n+1 


n+l 


n+l 


n+l 


p . = D u+E  v+F  p + G t+H 
Pn+!  Pn+1  n pn+1  n Pn+1  n Pn+1  n Pn+1 


t = D^_  u+E.  v+F.  p + G,  t+H 
n+1  tn+i  n tn+l  n fcn+l  n bn+l  bn+l 


(F7) 
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The  quantities  D,  E,  F,  G and  H are  obtained  at  all  the  grid 
points  by  the  following  procedure.  By  substituting  (F7)  into 
(F3) , (F4)  and  (F5) , to  eliminate  variables  with  subscripts  n+1, 

the  following  equations  are  obtained, 


bllun  + b12vn  * b13pn  + bHln  * dl  ‘ allun-l  ' a13Pn-l  <F8> 


b21un  + b22vn  + b23pn  + b24bn  * d2  ~ a21un-l  ' a23pn-l  ' a24  Vl 


and 


b,,u  + b,_v  + b,,p  + b,.t  = d, 
31  n 32  n 33*n  34  n 3 


(F9) 

(F10) 


where 


b. , = b. . + c. . D„  + c.,D 
1!  11  11  un+l  13  pn+1 


b.  t ~ b.  _ + c,  , E + c. 

12  12  11  un+l  13  pn+l 


b13  * b13  * CllFun+1  * CUPpn+1 


b!4  - b!4  * bUGu„+1  " b13=pn+1 


1 1 11  un+l  13  pn+l 


b21  ~ b21  + C21°un+1  + °23Dpn+1  + C24°t 


n+1 


28  b - ^ + c*i E + c ^ fjE  + c«*E. 

22  22  21  un+1  23  PR+1  24  tn+1 


b23  = b23  + C21Fun+1  + °23Fpn+1  + C24Ft 


n+1 
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b24  ■ b24  + =21Gun+1  + =23Gpn+1  + =24%^ 


d„  = d^  - c,,H  - _ - c,,H 

Jn+1 


2 2 21  u_  , , “23"p  _!.1  "24“tn+i 


b~,  = b,,  + C,,  D + C-.,D  + c„D  + c_.D, 

31  31  31  un+1  o2  vn+1  33  pR+1  34  tn+1 


b,_  = b,_  + C-.E  + c,_E  + c,,E  + c_. .E. 

32  32  31  un+1  32  vn+1  33  pn+1  34  tn+1 


b__  = b,_  + c,,F  + c,.F  + C-.-.F  + c,  ,Ft 

33  33  31  u_ , t 32  v_ , , 33  p„  . , 34  t 


n+1 


n+1 


n+1 


n+1 


b-,.  = b-,,  + c-.G  + c.0G  + c^G  + c^.G,. 
34  34  31  u_,,  32  v_ , , 33  p.,,  34  t 


n+1 


'n+1 


n+1 


n+1 


d_  = d,  - C.,,H  - - c,_H  - c-j  .H . (F 11) 

3 3 31  un+1  32  vn+]_  33  pn+1  34  tn+1 


By  solving  equations  (F6) , (F8)  , (F9)  and  (F10) , the  variables 

u , v , p and  t can  be  obtained  in  terms  of  u - , v , , p , 
n n ^n  n n-1  n-1  rn-l 

and  tn_j , and  in  the  following  form 


u = D v ,+E  V , + F p . + G t , 
n u n-1  u n-1  u n-1  u n-1 
n n n n 


v=Du  - + E v , + F p , + G t . 
n v n-1  v n-1  v n-1  v n-1 
n n n n 


p=D  U i + E v .,+Fp  , + G t 

rn  p n-1  - "-1  *'  ” 


£ n 


p n-1  p n-1  p n-1 

rn  rn  *n 


t = D.  u . + E,  v . + F,  p , + G.  t , 
n t n-1  t n-1  t mi-1  t n-1 
n n n n 


(F12) 
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where  the  coefficients  etc.,  are  given  by, 

n n 


Du  = allDll  + a21D21  + a41D41)/Det 
n 


Eun  * (a42D41>/Det 


Fu  " a13DU  + a23°2l  + a43D41)/,Det 
n 


Gu  = a14Dll  * a24°21  * a44D4l)/Dat 
n 


Hu  3 (dlDll  ' d2D21  + d3°31  ' d4D41)/Det 
n 


°v  = (all°12  * a21D22  - a41D42)/I)at 
n 


Ev  * a41D42,/Det 

n 


Ev  * (a13°12  - a23°22  - a43D42,/Det 
n 


Gv  * ,a14D12  ' a24D22  - a44D42,/Det 
n 


«v  - <-  dl512  +.d2D22  - d3D32  + d4D42)/Dat 
n 


Dp  = all°13  + a21D23  + a41D43)/Det 


Ep  = ta42D43>/Det 


F„  * ("  a13°13  + a23B23  + a43D43)/Det 
- n 
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Gp  - <’  a14513  + a24°23  + a44543)/Dat 

Hp  - iaiE13  - a2523  + a3S33  ' d4643l/Det 

Dt  - (all514  - a21624  - a41544,/Det 
n 

Et  = (‘  a42544)/Det 
n 

Ft  = U13514  - a23S24  * a43644,/Det 
n 

Gt  = a14°14  “ a24°24  " a44°44)/Det 

n 

Ht  = ai514  + 52524  " 53534  + d4544)/Det  {F13) 

n 

and  the  coefficients  6^,  D21»...  etc.,  are  given  by  the  following 
determinants , 


D11  = 

b22 

b23 

b24 

b32 

b33 

b34 

b42 

b43 

b44 

531  = 

b12 

b13 

b14 

b22 

b23 

b24 

b4  2 

b43 

b44 

' D21  = 

b12 

b13 

b14 

b32 

b33 

b34 

b42 

b43 

b44 

' °41  " 

b12 

S13 

b14 

b22 

b23 

b24 

b32 

b33 

b34 
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512  = 

b21 

b23 

£24 

' D22  = 

bll 

b13 

b14 

£31 

b33 

b34 

b31 

E33 

b34 

b41 

b4  3 

b44 

b41 

b43 

b44 

532  = 

bll 

£U 

b14 

' °42  " 

bll 

b13 

b14 

621 

£23 

b24 

b21 

b23 

b24 

b41 

b43 

b44 

b31 

b33 

b34 

613  b 

b21 

b22 

b24 

' 523  3 

bll 

b12 

b14 

£31 

£32 

£34 

b31 

b32 

£34 

b41 

b42 

b44 

b41 

b42 

b44 

633  = 

bll 

b12 

b14 

' D43  = 

£11 

b12 

b14 

b21 

b22 

b24 

b21 

b22 

b24 

b41 

b42 

b44 

b31 

b32 

b34 

514  “ 

b21 

b22 

b23 

' 524  = 

bll 

b12 

b13 

b31 

£32 

£33 

b31 

b32 

b33 

b41 

b42 

b43 

b41 

b42 

b43 

634  = 

bll 

b12 

b13 

' °44  = 

bn 

b12 

b13 

— 

b21 

b22 

b23 

E21 

b22 

b23 

b41 

b42 

b43 

b31 

b32 

b33 
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and 


Det  = 

B11 

b12 

513 

b14 

b21 

b22 

b23 

b24 

b31 

b32 

B33 

b34 

b41 

b42 

b4  3 

b . . 
44 

Equation  (E13)  gives  the  necessary  relations  for  the  coefficients 

D , D , D , etc.  Given  the  coefficients  at  n=N,  the 
u v p 
n n *ri 

quantities  at  all  grid  points  can  be  calculated  with  a sweep 
from  grid  point  n=N-l  to  grid  point  n=2.  The  values  at  n=N 
are  evaluated  using  equation  (E2)  which  gives. 


D =E  = F = G =0,  H =1 
“n  “n  ^ UN 


D = E = F = G =0,  H =1 


•N 


rN 


■N 


rN 


D.  = E.  = F.  = G,  = 0 , H = 1 


^ ^ *1* 


"N 


Dvn  _ " a41/b42  ' Evn  “ a42/b42 
Fvn  = " a43/b42  ' Gvn  = - a44/b42 


and 


H 


v. 


N 


(<34  “ b41  “ b43  PN  ’ b44  fcN)//b 


42 


(F14) 


Once  all  the  coefficients  of  the  recursion  formulas  are 
obtained,  the  solution  of  the  difference  equations  (32)  to  (35) 
are  obtained  with  a backward  sweep  from  n=l  to  n=N-2  using 
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equations  ( E 7 ) , provided  that  u^,  v^ , t^  and  p^  are  known  from 
(El)  and  given  by. 


and 


Pi  = £(d3  ’ 

C _ - H 

31  u2 

c32  Hv2  ' 

O H - 

j 3 P2 

- (b34 

+ c - . G 

31  u2 

+ =32  % 

+ G 

33  p2 

+ =34  S4,twl/ 

(b33 

+ c_,  F 

31  u2 

+ =32  Pv2 

+ CV1  Fr, 

33  p2 

* =34  Ft2> 

In  Appendix  G,  the  subroutine  DEQSO  of  the  computer  program 
provides  the  solution  of  the  present  finite  difference 
equations. 
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APPENDIX  G 
COMPUTER  PROGRAM 


The  following  computer  program,  written  in  Fortran  IV, 
was  used  to  solve  the  full  viscous  shock  layer  equations  as  a 
coupled  set.  The  input  quantities  are: 


Main  Program 


RMAC  Free  stream  Mach  number,  M . 

oo 

BO  Wall  to  stagnation  temperature  ratio,  tw/t  . 

REYIN  Free  stream  Reynolds  number,  Re^. 

TEMP  Free  stream  temperature,  T^  in  °R. 

ITH IN  0 for  thin  layer,  1 for  full  layer. 

IE  Number  of  mesh  points  in  the  n-direction. 

I END  Number  of  mesh  points  in  the  s-direction. 

IPRINT  Number  of  s steps  between  successive  printings,  s <_  1.0. 

IPRINl  Number  of  s steps  between  successive  printings,  s _>  1.0. 

DS  s step  size. 

DT  Time  step  size. 

GAM  Ratio  of  specific  heats,  y. 

SIGM  Prandtl  number,  a. 

XFACT  Convergence  criteria  for  local  iterations. 


Block  Data 

XNSH  Initial  shock  thickness,  assumed  constant  over  the 

body  (0.1072) . 
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Output  Quantities; 


MINF 

TINF 

TW/TO 

PR 

S 

X 

R 

NSH 

XSH 

RSH 

EPS 

NO  ITER 

NTIME 

CF 

HEAT 

STAN 

PWALL 

PW/PO 

uus 

ws 

PPS 

TTS 

RRS 

N/NSH 

CJ/USH 


Free  stream  Mach  number. 

Free  stream  temperature. 

Wall  to  stagnation  temperature  ratio,  t^t^. 
Prandtl  number,  <j. 

5,  surface  distance. 

Axial  distance  measured  from  the  body  nose. 

Body  radius  measured  from  axis. 

Shock  stand  off  distance  normal  to  body  surface. 
Shock  axial  distance  measured  from  body  nose. 
Shock  radius  measured  from  axis. 

Defined  as,  [p  (U  /C  ) /p  U a ] x/ 

co  p oo  OO 

Number  of  local  iterations. 

Number  of  time  iteration  cycles. 

* * * 2 

Skin  friction  coefficient,  2r  /p  U 

W w co 

* * * 3 

Wall  heat  transfer,  qw/pBUB  . 

Stanton  number,  q^/fH^-H^). 

* * * 2 

Pressure  at  the  wall,  p /p  U 

Pressure  at  the  wall/Nose  stagnation  pressure, 
u-component  of  velocity  behind  the  shock, 
v-component  of  velocity  behind  the  shock. 
Pressure  behind  the  shock. 

Temperature  behind  the  shock. 

Density  behind  the  shock, 
n coordinate. 

u;  normalized  tangential  velocity. 
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V 

P/PSH 

T/TSH 

R/RSH 

MACH 

PITO 

List  of 

SHVALS 

GEOM 


DEQSO 

DERIV 

MANISH 

BLOCK 

DATA 


Velocity  component  normal  to  the  body  surface. 
Normalized  pressure. 

Normalized  temperature. 

Normalized  density. 

Mach  number. 

Pitot  pressure. 

Subroutines 

Calculates  properties  behind  the  shock. 

Calculates  body  geometry  for  any  given  longitudinal 
location,  s.  This  subroutine  is  only  for  hyperboloid. 
The  asymptotic  angle  is  45a.  It  can  be  changed  by 
correcting  the  second  statement  (ANG  = . . . ) . 

Solves  coupled  set  of  equations  simultaneously. 
Calculates  the  derivatives  for  any  shock  shape. 
Calculates  the  new  shock  shape,  final  sweep. 

Ininital  shock  shape. 
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DIMENSION 

DIMENSION 

DIMENSION 

DIMENSION 

DIMENSION 


) 


1001 

1002 

1003 


900 


901 


COMMON/DEQS1/U1<111>  ^1<1U>»P1<111>-  rtdll  j-UC(  lll)^C<lll)* 

1 fc<  m ) , rc< m > ,uin( til > »vin< l n > ,pin< n i ) , riN< it l ) » uln<  n t > , 

2 VCN(lll),  F‘CN  (111)  jrTCN(lll)  ,U1NN(11  L ) ,TiNN(lit  > 

3 RKLll>*R2<lil)»RC(lll>  r.  X 

DIMENSION  VTSC(lll)  ,RVISC(  111  > ,RNSH(  111  > , RC3F  <111)  ,RCSFP<  111)? 

i PFAC< 111 ) ,XM< 111 ) »PITO< 1 11 ) 

us(iii)s>vs(iii>»ps<iii),rs<iii> 

F'O  ( 1 1 1 ) t PON  C 1 1 1 ) * PA  < 1 1 1 ) » PAN  (111) 

P 21(111)  »F‘21N  < 1 1 L ) , F'3'3(  1 1 1 ) , P33N ( 1 1 1 > 

00  < 1 1 1 > » VON (111),  VG (111), VGN  <111) » VOS (111) * 001 <111* 

OVSG  ( 2 ) , VSF'P(  2 ) 

COMMON/DEQSS/  U2< 11 1 ) » V2< 111 ) »P2< 111 ) »T2< 111 > »U2N< 111 > » W2N< 1 11 ) » 
i P”?N<  111 ) r T2N<  111 ) >U2NN<  1 1 1 ) »T2NN<  1 11 ) 

C0MM0N/BeQS/AllClll>.A13<lH>.BUCllt>.SI2<lll>.B13(lliWB14Clll  . 

1 C11(111>»C13<111)  , Dl<m>»A21(lil>*A23<lll>  »A..4  111>*B--1<|1 

2 ( 111)  , B23  (111)  ,B24(U1)  ,C2I  ( 111 ) ,C23  < 1 1 1 > * C-4  < 1 1 1 ) * l ’ 

3 B31(111)»B32<111) »B33<111) »B34(lll)»C3l<lll).C32(llt)»C33(lll>» 

4 034(111),  D3 (lil),A41(lll), A42( 111) ,A43(111) yA44(lll), D41 (111), 

5 B42(ill) , B43 (111) , B44 (111),  D4<111) 

~ CGMMON/MA INN/CRN I ,DS,DN< 111 ) , XN(111)*IM,I£ 

COMMON/SHCKG/XN5H  ( 202  ) , XNSP  < 202  ) , XNSF'P  ( 202  ) , VNSH  ( 202  ) , YNSP  ( .0^  ) , 

1 YNSPP  ( 202  ) » AD  (202) , AB(202)  > AC  (202) 

COMMON  /INSH/ 

1 

COMMDN/OUTSH/ 

1 

2 
3 

DIMENSION  YYYYC202),  XXXX(202> 

DATA  DFULL/ 'FULL- ' / » BTHIN/ ' THIN- ' / 

READC5, 1001 ) RMAC , DO ,REYIN, TEMP 
READ (5, 1002)  TTHIN, IE , IEND > IPRINT , IPRIN1 
READ < 5 , 1003)  DS,  DT  »GAM  , SIGM  , XFACT' 

FORMAT (4F10 . 3) 

FORMAT (514) 

FORMAT (4F7. 3, DIO.  3) 

CRNI  = 1 * ODO 
THIN=1THIN 

WRITE<6,?00)  y 

FORMAT  ( 1H1 , 45X,  'SHOCK  LAYER  PROGRAM  (U.  HOSNY  ) //3X, 

U)RITE(6,901 ) RMAC  , TEMP»REYIN  , BO  ?GAM,  SIGM  v 

FORMA T < 6X » ' MI NF ' , 6X , ' TINF ' , 6X , ' REYIN ' , 7X , ' TW/ 10  , 5X , GAM  ^7  .* 

t ' PR  ' ,/5X,F6.2»4X»F£.2»4X»F7.2»AX»F4.2,6X»F3.1»6X,F4,-) 


CONO 

t 

GAM 

9 

S » 

UPSH 

9 

XNS 

EPS 

r 

RMAC 

9 

TF’SH  , 

VI  SCO 

PPS 

f 

RRS 

9 

TTS  * 

UUS1 

9 

VOS 

PPS1 

9 

RRS1 

y 

TT51  , 

UUS2 

9 

UVS1 

PSP 

f 

RRS2 

9 

TSP  , 

USF 

r 

VVS2 

PPS2 

9 

RSP 

9 

TTS2  , 

UUS 

4 

VSP 

INPUT  DATA'/) 


CNT=BTH1N 

IF  < 1THIN ►£□.  !>  C NT'S FULL 
UIRITE<  6, 902  > CNT  , IE  , IEND  » DS,  DT 
902  FORMAT ( /3XA5 , ' SHOCK  LAYER ' » 5X , ' NO . 

1 'NO.  OF  STEPS  IN  3 ='»I4,4X, 


OF  STEPS  IN  N = 
'S  STEP  SIZE  = ' 


,14, 4X , 
F5.3,4X, 
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2 ' T STEP  SIZE  = ' y F6 , 1 > 

DO  LO  ! -i > TEND 

10  xNsr < t ) -0 * ono 

IM-1E-1 
Y1M- CM 

nr- 1 ; jo  o / y i h 

xn  <1)^=0.  ono 

BG  J.5  N=1*IM 
TiN  < N > — DY 

15  XN<NtL>-XN(N)+J)N<N> 

YNSH ( 1 ) =0 . 0 

CALL  fiERIV ( HGs>  lENDy  t ) 

IENm«IENn+l 

WR  f TE ( 6 » 903  ) ( XNSH  < I ) » 1 = 1 * TEND  L ) 

WR ITE <6,904 > ( YNSII(1>*  [-tr  IENDL) 

903  FORMA  T < /3X  > ' INITIAL  SHOCK  SHAPE  ' , ■'  / 5X  »'  SHOCK  THICKNESS'/ 

1 30  < 3X  ■>  1 OF  1 0 » 4/  > > 

904  FORMAT <//3X, 'SHOCK  KAIUJ  IS ' /30 t 3X , 10F10 * A/ ) ) 

NTXME=0 

NT I ME 1-0 
TIME=0>0 
20  NT  IME-NT1ME1 .1 
TiME=rihEfnr 
RSH=0 » 0 
UPSH=0  , 0 
TPSH-0 , 0 
oiscd=o*opo 

CONQ=O.ODO 
XNS=XNSH<1  ) 

XNS L=XNS 

CMS=(XNB1 lXNS>/2. 

DS2*ns/2.orio 
ck=i .ono 
csf=o.odo 

STF=  L * ODO 

PIIIC=DARC0S(0.01!0) 

rs-o  ♦ ono 

RS2-O.ODO 
XB=Q » 0 
cnF=o . o 
CDP=0 « 0 
CtiP  L -0,0 
crip  2^0.0 
CDF 1-0,0 
CDF 2 = 0.0 


CDPD-0 .0 
CDFD-0 . 0 

PO  T P - < ( GAM  +1  ,0>*RMAC*RHAC/2.0}**<UAM/',  GAM-  1,0)  )/<OAM  HcMACjKKMAC# 

1 <2.0*l)AM*RJlAC*RhAC/  ( GAM!  L . 0)  -(GAM-- 1 . 0 ) / (GAM  1 1 .0  ) ) 

2 ( 1 , 0/ ( GAM-1 , 0 > ) ) 
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TW=BO*< 1 .0/<  CGAH-l .0)*RMAC*RMAC) fO.50) 

TB-TUI#(  (GAM-1  .0)  *RMAC#RMAC*TEMI- > 

CGNP=198.A/<  (GAM-1.  ) *RMAC*RMAC* TEMP > i 

VISRA=  C 1 .OFCONP}#<  1 . 0/  ( (GAM-1  ► )*RMAC*RMAC)  ) **1  . S/  < i . / ( (GAM-1  . > * 

1 RMACflRMAC  ) 1C0NP) 

EPS-1. 0/»SCIRnVISRrt*REY  IN) 

CALL  SHOALS  < 1 . ODO,  0.000* 1 .000*  0.000*11  SO*  WSOjtJUSO*  PPSO.  1 ) 

TTS=TTSO 

DO  100  N=1*IE 

RNSH(N)=CNS/( 1 ,+CK*CN3*XN(N) > 

RCSF ( N ) -CNS/ ( 1 . FCK*CNS*XN ( N ) ) 

RCSFP ( N > = 1 . 0+CNS#XN ( N ) 

U1 (N)=XN(N) 

U2 ( N ) — XN ( N ) 

U1N ( N > =1 .0 
U2N(N>=1 .0 
U1NN ( N > =0 . 0 
UC ( N ) =XN  < N ) 

UCN(N>-1 .0 
US  < N ) =0 . 0 
OS  < N ) =0 • O 
F'S  ( N ) =0 . 0 
TS  < N ) =0 . 0 
01 ( N ) =XN ( N ) 

02(N>=XN (N) 

OC(N) =XN(N) 

Ti<N>  = 1.0-U.0-XN<N> >*<i.0-TW/TTS0) 

T2  < N > =T1 ( N) 

T1N(N)=1 .O-TU/TTSO 
T2N ( N ) =T IN ( N > 

T1 NN ( N ) —0 . 0 
TCCN)— T1(N> 

TCN(N)=T 1N<N) 

UISe(N)=<TTS-t-CONP)*TC(N>**l . 5/ < TTS*TC < N > +CONP ) 

RVI  SC ( N ) = ( TTS*TC ( N > + 3 . O&CDNP )/<2»0#TC(N)*<  TTSfcTC (N)4CQNP) ) fcTCN ( N ) 

PI (N>=1 .0 

P2 (N ) =1 .0 

PC ( N > =1 . 0 

F'A(N)  = 1 .0 

PO ( N > = 1 ♦ 0 

PON (N) =0.0 

RC(N)=PC(N)/TC(N) 

FCN ( N ) =0 . 0 
PI N ( N ) =0 . O 
PAN<  N)-0. 0 
P33(N) =O.OKO 
P21 < N ) =1 . QDO 
100  P2N ( N )=0 . 0 

OISCO=( 1 .0+CONP)*TTS**1 . 5/ ( TTS+CONP ) 

CONO=VISCG/SIGM 
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DO  5000  T~  L , I END 
YT-I 

,ODO)+DS 

CALL  GROM  < 3 , DS2 . PS .2 , CK2  , 331-2 .3  [1-2,  XB2 ) 

IFt  I.Efl.l)  CK1=CI\2 
CKT-(CK2-CK1  ) / DS 
PH  IP*  -Cl\ 

F'HI  = DARCOSICSFS) 

XNSF‘M~  < XNSF'  ( I > PXNSF't  IM  1.W2.0D0 
CNSN=XNSH(I> 

XNGN- ( XNGH < I > FXNSH< I FI ) )/2.0 
[Fd.EQ.D  CNSN-XNSN 

ALP  - PHI  + riArAN(XNSPM/<l,0H:i\2*XNSN>) 

ALh'C -PH IC+DAT AN < XNSP < I )/ < l*0+CK*CNSN)  > 

S P = DSIN(ALP) 

CP  * DCOS(ALP) 

SPH-SP*SlF21-C|,,'*CSIr2 
CPB~CP#0  JIF2-BP#CSF2 
IFCS.Lr+O.OOQl > GO  TO  120 
RSCP*<RS2-RS1  )/IiO 
CSFF’=  ( CSF2--CRF1 ) /D3 
120  CONTINUE 
NITT-0 

2000  N I n =NITTH 

CALL  SHOALS  ( SP , CP,  SPB,  CPU,  TTSH , VRSH,  URGH,  PPSH,  2.1 

usriausp 

RSDaRSP 

GAMQ=GAH*GAM 
GAMH-GAH-1 .0 
GAMB=GAh+l » 0 

GAMT-GAH*GAMB/CAMM-2.0*tjAM*GAMM/GAHB 
CPHC=DCOS(PH  tC> 

SPC=DS1‘N  ( ALPC  > 

CPC-DCGS ( ALPC ) 

S2PC=2 . 0*SF'C*CPC 

SPCQ=SPC*SF‘C 

SPCT=SPCO#SPC 

S2PF‘H -USIN  < 2 , 0* ALPC -F’FII C) 

SPPH=DSIN<ALPC--PH  CO 
CPF’H=DCOS  ( ALPC—F’HIC ) 

C2PF'H=DC0S  ( 2 . 0*ALPC  -PHIC  ) 

RMACQaRMAC*RMAC 

AR=2  ► 0*0AM*RMACO*SPCQ  -GAMM 

AK1P=2 .0*GAMtl*RMACG*S2PC/ < AR*GAMB > 

1 F4 . 0*GAM*CPC/  < GAMB*RhACG*AR*SPCT> 

-4.0*GAM*GAMtnRMACG*RMACa*S2PC*SPC()/(ARYAK'*l!AMEO 
- 2 > 0#GAM*HMACQ*S2FC*0AHT/  ( AR*AR  > 
+4.O*GAMO#y2PC/(0AMD*AR*Ak!KSPCQ) 

AKla-(  1 ,0-<GAhH/2.0)H<Trs/PPS>)*a2F,PH+GrC4SPPH*BAMM*AKLP/fiAh 


85 


AEOC-TR-79-25 


AK2=CPC*SHPH-  ( GAMM/GAM ) * C TTS/PPS ) *SPC#CF'PH 
AK3=2  . 0#S2F*C/GAMEt 

AK4=2.O*GAM*S2PC/<GAMB*GAMB>-M.0#erC/<GAMB*GAMF«*RMACn#RMACG*SF'Cr> 
AK5=  < 1 . 0-  ( GAMM/GAM ) * < TTS/PPS  ) > #C2PPH-SPC#CPPH*AK1  P*GAMM/GAM 
AK6=-CPC#CPPH-  SPC*SPF’H#  ( GAMM/GAM ) * < TTS/PPS ) 

ALPF'=CXNSPP<  t >-(CNS-CNSN)#2.0/<3.0*DT>  >/<i.0+CNS  )-l  .0 

USF'=0.0 

PSP-0 ♦ 0 

TSF'=0 . 0 

RSF-0.0 

IF(I.EO.l)  GO  TO  123 

ALF‘F'=  < CPPH#CPF'H*  ( YNSPP  ( I ) - ( RSH-YNSH  i I ) > *2 . O/DT  ) -YNSP  C 1 ) * < CK*2 . 0* 

1 SPPH*CPPH+CNS  *CKP*CPPH*CPF'H/C1.0FCK*CNS  )))/<< 1 . 0+CK*CNS  >* 

2 CPHC) 

USP=AK5*ALPF+AK6*FHIP 

PSF'=AK3*ALPF' 

TSF-AK4*ALF‘P 

RSP=  ( GAM/GAMM  ) * ( PSP#TTS-TSP*PF'S ) / < TTS*TTS  > 

123  CONTINUE 

USP=AK  1 # ALPP-f  AK2*PH  TP 
C 

IF(I.EQ.l)  WM=WS 
IF(I.EQ.l)  VSn=VSP 
I F < I . NE . 1 ) VUM  = 1»<>00 
IFCI.NE.l)  VSn=0,0D0 
DO  110  N=1 » IE 
IFCS.GE. 0.0001)  GO  TO  108 

PFAC  ( N > =4  ♦ 0D0*  CPA  < N ) + ( PF  S2/PPS0-2  . ODO ) #F'0  ( N ) ) / < UUS2*US ) 

1 -XNSP(2)*XN<N)*F'0N(N)/(2.0n0*UUS2*CNS) 

GO  TO  109 

108  PFAC<N)  = (PS(N  >-XNSP(I)*XN<N)#F  CN(N)/CNS+PSP*PC<N>/PPS>/UUS 

109  CONTINUE 

RNSH(N)=CNS/< 1 . +CK*CNS*XNCN> ) 

UISC(N)  = < TTS+CONP  >#TC(N)4<#1  . 5/ < TTS*TC  ( N ) FCONF' ) 

RUISC  < N ) = < TTS*TC  < N ) +3 . 0*C0NP ) / ( 2 . 0#TC  < N ) * < fTS*TC  C N ) FCONP ) > *TCN ( N > 

110  CONTINUE 

VISCO=(  1 .0+C0NF')*TTS**l . 5/ < TTS FCONP > 

CONOCO I SCO/ S I GM 

REFAC— RRSfcOUMfcCNS/  C EPS*EPS*U  ISCO  ) 

IFCS.LT.O. 0001 ) GO  TO  170 
XNSF'T=(  XN5-XNS1  )/DS 
DO  165  N=1 r IE 

RCSFCN)=CSF*CNS/(RS+CNS*XNCN)*CSF) 

RCSFPCN)=RSCP+XN(NmXNSPT*CSF  + CNS*CSFP> 

165  CONTINUE 
170  CONTINUE 

C S MOMENTUM  COEFFICIENTS 
DO  400  N=2*TM 

Ai=-REFAC*RNSH(N)*U(JS*RC(N)*UC(N>/<  VISC ( N ) *UUM ) 

B1=REFAC*(RNSFKN)*UUS*XNSFC I )*XN  ( N)  *RC  ( N ) *UC  ( N ) / C UVM*CNS*UTSC ( N ) ) 
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1  -RC(N>#UC(N>/UISC{N>  > FRUISCCN)  FCKrtRNRH  < N > + RCSF  < N ) 

Cl  = H^NSH(N)*<REFAC*RC<N)*U2.0*U5F*UC(N>/(  WMAUTSOtN)  )+CK*UC<N) 

1 /UISC(N)-UUS*XNSP<  I )!KXN<N)*UCN(N)/<WH*nNS*VISC(N)  HUUG*US<N>/ 

2 <UUM*UISC<N) ) )+CK*<RVTSC<N)  tCKY-RNSH  C N ) +RCSF < N ) ) ) 

C2=-REFAC*RC  ( N ) * ( UCN  ( N ) /V ISC  < N ) +C!<  *RNSH  ( N ) *UC  ( N ) /V  ISC  ( N ) > 
C4a-REFAC*RC  ( N ) * ( RNSH  < N > * ( UUSflXNSP  ( i ) *XN  < N ) HJCN  < N - / ( UUM*CNS* 

1 uisc(N)  )-usp*uc(N)/(yyh*y.rsc(N)  >-uus*us<N)/< wm#visc<n>)-ci\* 

2 UC(N>/VISCCN>  )*UC(N>  -UC  ( N ) *UCN  ( N >/(J  ISC  < N ) ) /TC  CN  > 

D=REFAC*  < RNSH  C N ) * < ( -UUS*XNSP < t ) *XN  ( N ) *(JCN < N ) / < UWWJNS  ) +UQPSOJC  < N ) 

1 /WM+UUS*USCN)/WM>*RC(N>*UC<N>/UISCIN>+CI\*RC(N>*UC  (N)*UC<N>/ 

2 VTSC(N)  >+RC(N)#VC<N>*UCN(N)/UISC<N)  ) 

IF  < S . GE . 0*0001 ) 00  TO  330 

A3-0.O 

B3=0.0 

C3=-C4/RC(N> 

D=Ei-REFAC#RNSH  ( N ) *PPS*PFAC  ( N ) / ( RRS#WM*U  ISC  ( N > > 

GO  TO  331 
330  CONTINUE 

A3=  -REFAC*RNSH { N ) #FF'S/  ( RRS*IJUS*UUM*VISC  ( N ) ) 

B3=ftEFAC*RNSH ( N ) *PPS*XNSF‘ ( I ) *XN ( N ) / ( RRS-WJUG*  WM#CNS*UISC ( N > ) 
C3--C4/RC < N ) -REFAC*RNSH < N > *PS P/ ( RRS*WM*UIJS*UISC ( N ) ) 

381  CONTINUE 

EiNl— DN<  N)+DN<  N— 1 ) 

DN2»BN<N> ~BN(N-1 ) 

All <N>=(2.0~B1*DN<N) )/BN < N-l )/DNl 
A13<N>=-<B3*rtN(N)  )/DN(N-l  )/HNl 
Btl<N)«Cl+(-2.0+Bl*DN2>/DN(N)/nN(N-l )+Al/<BS*CRNI ) 

B12 ( N )~C2 

B13(N>=  C3  + B3*EiN2/DN(N>/DN(N~l > +A3/ < BS*CRNT  ) 

B14 ( N ) =C4 

Cll <N>“<2,0+Bi*DN<N-l> )/DN<N>  /UNI 
C 1 3 < N ) -B3*DN ( N- 1 > /DN ( N > /ON 1 

D1(N)=CA1#U1 CN>+A3#P1 (N> )/( BStfCRNI ) -H/CRNT- ( ( L . 0 "CRN I ) /CRNI ) # 

1 CU1NN<N)+BI*U1N<N)+B3*P1N<N>FC1*UI (N) FC2*UL(N) FC3*P1 (N)+ 

2 C4#T 1 ( N > ) 

400  CONTINUE 

C ENERGY  COEFFICIENTS 
JOO  420  N*2  t IM 

A3=REFAC*RNSH<N)*UUS#PPS*UISC0*UC<N>/<WM*RRS*rTS*C0NG*Ul>SC(N)  ) 
A4'=-REFAC*RNSH(N>*UUS*UISCU*RC(N>*UC<N)/<WM*C0N0*UiSC<N>  > 

B1  =2  * 0#UUS*UUS*Vi: SCO# I UCN  < N > -CK*RNGH  C N ) #UC ( N > ) / < TTS*CONO ) 
B3==REFAC*VI  SCO*F'F'S*  < -UUS*RNSH  ( N ) *XNSP  ( I ) *XN  ( N ) *UC  < N > / < WH*CN S.i  ) 

1 +UC(N) >/ICONO*TTS#RRS*UISC(N) ) 

B4-REFAC#  < UUSfcXNSP  ( J ) *RNSH  ( N ) *V  [BCO*XN  < N ) #KC  ( N ) *UC  < N > / < WM*CNS* 

1 CONO*UTSC(N> >-VISCO*RC<N>*UC<N)/(CGND*','ISC<N) > ) FRUISC(N)+ 

2 CK*RNSH<N> FRCSF(N) 

Cl-RNSI-K  N)*UISCO#(REFAC*  ( RC IN)* ( 1JUSXXNSF’ ( l )#XN  < N>*TCN< N > / < UUM*CNS  ) 

1 -UUS#TSP#TC<N)/(UVM3KTTS)  -UUS3KrS(N)/WM)/VTSC<N)-»UUr>#lJUS;HPr,S'* 

2 PFAC<N)/<  WMXtRRS*  TTS#VISC<N)  ) ) -2 , 0#CN*UUS»UUS*<  LJCN  ( N ) -CKt 

3 RNSH<N)*UC<N> )/rrS>/CQNO 
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C2=REFAC*U  TSCO*  C -PC  < N ) * FCN  C N > / FC  < N ) FPPSYPCN  ( N ) / < RKS*  F FS  '>  ) / 

1 lCONO*UISC(N> ) 

C3-=REFAC*y  l sen# < RNSH  < N > *UC  < N ) * C < UUG*XNG P < I > #XN  ( N ) f CN  l N ) / < WM*LNS» ) 

1 -uug*tsp* re  ( n > / < whi  r rs  > -uus* rs < n ) zoom ) / re  < n ) hjus*psp/ < yyh*RKS 

2 #rTS>>-VC'N>*TCN<N>/TC<N>)/lCONO*VISC<N>> 

C4=-REFAC#U1SC0*  ( RNSH (N  )*h:C<N>*UC(N)*  UHJSttXNSP  ( I > ( N ) * FCN  < N ) / 

1 (yVH#CNS)-UUG#TS<N)/WM  J-RC(N)  *OC(N  ) * TCNtN) ) / <.  C0N0-*  TC  ( N ) 

2 *yiSC(N)) 

0  =-REFAC*U  l SCO*  i.  RNSH  ( N ) * ( RC  ( N ) *UC  < M ) # < UUSfcXNGP  < X ) *XN  < N ) * TON  t N > / 

( UVM*CNS  > -UUS*TSP*  TC  < N ) / ( VUM*T  FS ) -UUS*7  S ( N ) /UUM ) HJUG*PPS*UC  C N > * 
UUS*PFAC  ( N ) / ( uyM*RRS*TF  S ) ) -RC  ( N > *UC t N ) * TON  ( N ) +PPS*YC  < N ) *PCN  t N ) / 

< RRssrrs)  )/(cono*oisc<nj  > -uug*uus*u  xsco* ( ijcn<  n > -lk*rnsh<n ) * 

UC  ( N ) ) * (.  UCN  ( N ) -CK*RNSH  ( N ) *UC  < N ) ) / ( TTSfcCONG  ) 

DN1  = DN  < N > +DN<  N -1 ) 
riN2=EiM(N)-DNCN-i  ) 

A21 <N)=- (B1*DN<N) )/DN<N~l >/DNl 
A23<  IM)=-<B3*DN(N)  ) /DN < N -1 >/DNi 
A24(N>=<2.  0-B4*DN ( N ) )/DN (N-l )/»NL 
B21<N>=Ct+Bl*DN2/DNCN)/DN<N-i ) 

B22  ( N > = C2 

B23 ( N ) ~C3  FB3*DN2/EiN<N)/DN<N-l  ) FA3/(BS*CRN  l ) 

32 A < N ) -A4/  < DS*CRNI ) +B4*DN2/DN  < N ) /DN  ( N-  L ) -2 . 0/BN  < N ) /DN  ( N - L ) +C4 
C21  (N)=Bl*riN  (N~t  >/HN(N)/DNl 
C23<N)=B3*DN CN-1 >/DN<N)/HNL 
C24(N)=<2.0+B4*DN(N-1 ) >/DN<N)/DNl 

Ei2  < N)  = ( A3*P1  ( N ) +A4*T1  (N  ) ) / ( DStCRNI  ) -D/CRNI  -<  ( L » 0 -CRN  I > /CRN  [ ) * 

1 (TiNN  ( N)  +EU*U1N  (N  ) FB3*F’1N<  N)  TB4*T  LN  (N>  FCltUl  <N)  FC2*U1  < N>  + 

2 C3*F'l  X N > fC4*T  1 < N > ) 

420  CONTINUE 

C CONTINUITY  COEFFICIENTS 
DO  440  N=1 f IN 
UCH=(UC(N)+UC(N+l ) >/2»0 

och=  < oc  < n > +yc  c n+  j.  > ) /2  ► o 

F'CM=  < PC <N)  FF'C<N+1)  )/2.0 
TCM=<  TC(N) fTC( Nfl ) )/2.0 
XNM= ( XN  < N ) FXN  < N+ 1 ) )/2,0 
RNSHM=  < RNSH  C N > +RNGH < N+l ) )/2.0 
RCSFM- < RCSF ( N ) +RCSF ( N + i > )/2.0 
RCSFF'M=  C RCSFF’  ( N ) +RCSFP  ( N+ 1 ) )/2.0 
IJSM=  ( US  ( N ) +US  ( N F 1 > ) /2  * 0 
PSM=<PS(N)+RS(N+1 ) )/2.0 
TSM=(TS<N>+TS<N+1) >/2.0 
UCNM= ( UC  C N F 1 > ~UC ( N ) ) /BN ( N ) 

VCNM=< yC(N+l)-UC(N) )/DN(N) 

PCNM=<  PC(NF1  )-F'C(N)  )/DN(N) 
rCNH-<  TCCN+l ) -rC<N> )/DN<N) 

IF  (S.GE *0.0001)  GO  F0  430 

Al-0.0 

A3=0  * 0 

A4=0 ♦ 0 


88 


r-j  m -i  c j n ^ n n 


A EDC-TR-79-25 


01=0*0 

B2“RRS#VVM*F'CM#TCM 
B3=RRG*WM*VCM*TCM 
B 4 =-RR  SX V VM* VC  K*PC  M 
CD4.0*RNSHM*RRS#UUS2*PCM*rCM/DS 

C2=RRS*VVM*(2.0*RN!3HM*rCM#TCMf  IChXPCNh  -PCMXTCNMj 
C3-RRS*<2*0#RNSHM#TCM*(2.0*UUS2*UCM/DSFVVM*VCMHVVM*<  rCM*VCNi1-- 
1 VCM# TCNM) ) 

C4=RRS*  < 2 * 0*RNBHM#PCM#  ( 2 . 0*UU52#UCM/DS4  V VMXVCM  ) + VVh#  < PCM+VCNM 
1 + VCM*F'CNM>> 

D~RPS* < -4  >0#RNSHM#PCM#TCM* < 2 . 0*UUS2*UCM/DS FVVM*VCM 1 +2  * 0*0 VM* ( VCM* 

1 PCM*TCNM-VCM#TCM#F‘CNM-PCM*TCM*VCNM> ) 

GO  TO  431 

430  CONTINUE 

A1=RNSHM*KRS*UUS*TCM*PCM 
A3=RNSHM*RRS*UUS*UCM# TCM 
A4=-A3*PCM/TCM 

B1=-RNSHM#XNSP(I>*RRS*UUS*XNM*PCM*TCM/CNS 

H2=RRS#VVM#PCM*TCM 

B3=TCM*  < RRSfcVVMXVCM-RNSHMXRF^UIJStfXNSP  ( [ ) *XNM#UCM/CNS ) 
B4=-B3*PCM/TCh 

C1=RRS*UUS*RNSHM*  < -XNSP ( I > *XNM*  < TCM*PCNM~FCM*TCNM ) /CNS  f ( T'CMXP SM- 
J.  PCM#  TSM ) +TCMXPCM# < U6D/UUS+RSD/RKS FRCSFM* < RCSFPM/ ( CGFXCNS > - 

2 XNMXXNSF'  ( I ) /CNS  ) + < XNSPT-XNSP  ( I ) >/CNS)  ) 

C2=RRS*VVM*( TCM*PCNM-PCM*TCNM+TCM#PCM*CCK*RNSHM+RCSFM> > 

C3=RRS*VVM* C TCM*VCNM+ <CK*RNSHM+RCSFM ) *VCM*TCM- VCM* TCNM ) +RRS#UUS* 

1 RNSHM*  C TCM#  ( USM-XNM#XNSP  ( DXUCNM/CNS ) -UCMXTSM+  TCM*UCM*  ( UBD/UliS 

+RSD/RRS+RCSFM*  ( RCSFPM/  < CSF#CNS  > -XNh*XNSP  < I ) /CNS  ) + < XNSPT- 
XNSP  ( I ) >/CNS)+XNM*XNSP< I )#UCM#TCNM/CNS) 
C4=RRS*VVM*<PCM*VCNM+<CK*RNSHM+ftCSFM)*VCM*PCM-»  VCM*PCNM) fRRSXlJUS* 
RNSHM#  < PCM# < USM— XNM#XNSP < I > XUCNM/CNS > +UCM*PSM  t PCMXUCM* < USD/UUS 
fRSD/RRS+RCSFh# ( RCSFPM/ <CSF*CNS ) -XNMAXNBP ( I ) /CMS ) + t XNSPT  - 
XNSPC I>  )/CNS>-XNM#XNSF'<  I > *UCH*PCNM/CNS > 

0= -2  * 0*RRS*VVM* < PCM*TCM*UCNM+<CK*RNBHM+RCSFM ) *VCM#PCM* rCM+VCM* 

< TCM#PCNM-F’CM#TCNH ) ) -2 . 0*RRS#liUS#RNSHM*  < F'CM#  TCM*  ( USM-XNM* 

XNSP < I > #UCNM/CNS ) HJCM# ( TCM*PSM-PCM#  TSM ) +PCM#  TCM#UCM# < USD/ IJUS f 
RSD/RRS+RCSFM*  ( RCSFPM/ ( CSF*CNS > -XNM*XNSF‘  ( I ) /CNS)  + < XNSP  f- 
XNSP  CD)  /CNS  ) — XNM#XNSF‘ < I ) *UCM#  < TCM*PCNM~PCM*  TCNM ) /CNS  ) 

431  CONTINUE 

A41  (N+l ) -AH-CFCNI*  (DS#Cl-2  »0#B1#DS/DN  < N ) ) 

A42  ( N+l ) = , CRNI*(DB*C2-2*0*B2*DS/DN<N>  > 

A43  C MU  ) =A3  J-CRNI*  ( DSXC3-2 . 0#B3#DS/DN  ( N ) ) 

A44CN+1  >=A4+CRNI*(»G#C4™2.0#B4#DS/DN(N>  > 

041 (N+l )=A1 f CRN I#  C DS*C  t +2  * 0#B 1 XDS/DM C N ) ) 

B42CN+ I )=  CRNI#(DS#C2+2*0#B2#DS/DN(N) ) 

B43  CM+l ) -A3+CRNI * ( DS#C3+2  * 0#B3#DS/DN  C N ) ) 

B44CNH  >=A4fCRNI#<riS#C4  F2  * 0#B4*DS/DN  < N ) ) 

D4  < N+l ) --2  * O*0S#D+  ( A 1 — < 1 ► O-CRN I ) #C1#DS)  # ( U1  ( N M ) fill  (N)  )-<l  .O-CRN  L ) 

1 #C2#DS#  ( VI  ( N+l ) +V1  ( N ) ) + < A3-C 1 .0-CRNI  >#03#DS)X<  PI  (Nf.L  )+PUN)  ) 

2 KA4--C 1 . O-CRN  1‘ )*C4#HS)*<  fiCNf  L)  FT1  < N ) ) -2  * 0#  ( 1 ♦ O-CRNI ) # 
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3 < DS/DN(N)  )*(Hl*!Ui(NH)-Ul<N)  ) + B2*  < 0 l ( Nil ) -VI  (N3  3 4-B3* 

4 <P  L (N  + l >-F'l  (N>  )+B4*<  T1 (N  + l ) - T1 ( N 3 ) ) 

440  CONTINUE 

C N MOMENTUM  COEFF  tCIENTB 
DO  460  N=2  ? IE 
PCH-(PC(N>  1PC  ( N-.L  ) 3/2.0 
UCM=(UC(N) 1UC(N-1 3 3/2.0 
UCM— < OC  < N 3 !0C(N-l>)/2.0 
TCM-(TC(N3+TC(  N-.1  ) )/2.0 
XNM=(XN(N)+XN<N-1 ) 3/2.0 
RNSHM- C RNSH ( N ) 1RNSH < N-l 3 3 /2 . 0 
PCNM  - < PC(  N 3 -PC  ( ft-l >> /UN (N-l 3 
0CNM=<UC(N3-0C(N-1 ) 3 /BN < N-l ) 

0SM=(0S<N)+US(N-1 > )/2.0 

OUO=USM-XNSP< I >*XNH*OCNM/CNS  10Sn*UCM/00M 

C1=RNSHM*RRS*UUS#<-2.0*CK*UUS*PCM*UCM+T(+IN*OOM#OUU*PCM ) 
A2==RNSHM*RRS*00M*UUS*THIN*PCh*tJCM 

82-RRS*00M*  rHIN#FCM*  ( -RNSHMKXNSP  ( I 3 *XNM  *UlJS*UCM/CNS  + UOM*U(.:M  ) 
C2=RRS*00M*THIN*PCM*  ( RNSHM*USD*UUS*UCM/OUM  100M*UCNM  3 
g3_~ppc;jj<  TCH 

C3--RRSWON*<RNSHM#(THIN*UUS*OOU*UCM-CK*UUSWUS*UCM*IJCM/OOM3  + 

1  THINtfUUM*OCM*VCNM 3 
C4-PPS*PCNM 

H=  -PPS*TCM*FCNM+2 . 0*RRS*00M*  < RNSHM*UUS*PCM*UCM*  ( CK#UIJS*UCh/OUM- 
1 TH  IN'^OOO  3 -THI  N#OUM#PCM4tOCMiJiVCNM  3 
ECU  (N-l )~CRNI*DS*C1 

B32(N-1 )=A2+CRN [*<nS#C2-2.0*H2*DS/DN( N-t ) 3 
B33  ( N-l  )-CRNI*(HS*C3  -2.0*B3#DS/DNCN-1 3 ) 

B34(N-1 >=CRNl*nS*C4 
C31  CN-1>=CKNT*DS*C1 

C32 ( N -1 ) -A2+CRNI* ( DS#C2 12 , 0*E<2*D5/DN < N- 1 3 ) 

C33(N-1  )=CRNI*(  0S*C312.0*B3*DS/BN(N--1 ) 3 
C34(N-13=CRNl*riS*C4 

It3<  N-l  >--2.0*DS*DK  A2-<  l .G-CRNI  >*IiS*C2 ) # ( VI  (N)  10  L < N~1  ))-<!. O-CRNI 3 

1 *DS*(C1*(U1 <N 3+U1 (N-l 3 3+C3*(Pl(N3 1P1 (N-l ) )+C4*( T1(N)1T1<N-13)) 

2 -2.0*( 1 .O-CRNI  3#nS#(B2#  (01  (N3-V1  (N-l  3 )/DN(N-  l ) + B3'K 

3 ( PI  < N ) -PI  < N-l ) 3/LlN  (N-l ) 3 
460  CONTINUE 

T2 ( 1 )=TW/TTS/CRNI~< 1 .O-CRNI 3 *T1 ( 1 3/CRNX 
02IE=02(IE> 

CALL  DEOSO  ( TFACT»I3 
DO  802  N=1 » IE 
R2(N)=P2(N)/T2(N) 

802  CONTINUE 

IF ( Ni  TT . GT . 1 3 XNOO=XNSO 
XNSO=XNS 

IFtS.LE. 0.00013  GO  TO  217 
IF(NITT.EO.l')  XNS=1.01*XNS 
IF (N ITT  .Eli  ► 1 3 GO  TO  711 

XNS^XNS  1 ( XNQD-XNS ) * < U0S2-02  < IE  3 3 / < 02IE-02  ( IE  ) 3 
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GO  TU  7 11 
717  CONTINUE 
AA-0>0 
Bit' ••=0.0 

DO  700  N*»2>LE 

AA  =AA  +DN(N~l)K(R2<N-l)*U2<N“l>lFi2<N)*U2(N))/2* 

700  BB  “BB  HJN(N-1  > * < R2 l N -1 ) *U2 (N-l > *XN < N-l > 

1-  +R2(N)*U2( N)*XN<N) )/2. 

IF  ( S * GE  « 0 * 0001 ) GO  TO  705 
A£A=8.*BB  *RRS2»UUS2*CSF2/D3~DS 

B [ B-4  . *AA  #RRQ2*UUS2*Rn2/IiS-DS 

cic*-ns 

ROr=BIB*BIB-A,[A#CIC 

XNS  « (-BIB  1 DGQR'T  ( ROT 5 ) / AIA 

CG2  =*XNS*(RS2#AA  +XNS*CSF2*BB  >*RRS2*1JIJS2 

GO  TO  711 

705  AtA-DB  #CSF2#RRS2*UUS2 
B [B-AA  #RS2*RRU2*UUS2/2 . 

CIC--C01  +<RSfCNS#CSF)#<  (1  »+CK*CNS>#RRS*WS~XNSF,(  [ )#RR8#UUS>#DS 
ROr=BIB*BIB-AIA*CIC 
XNS  = (-BIB  + riSQRT(ROD)  / AIA 
C02  =XNS*(RS2*AA  +XNS*C3F2*BB  >*RRS2*IJUS2 

711  IF  ( S ♦ GE  » 0 » 0001 ) GO  TO  715 
XNS J -XNS 

715  CNS-CXNSFXNSt ) /2.0D0 
RSH1-=RS2+XNS#CSF2 
RSH  -RE  TCNS*CSF 
XSH=XB-CNS*SIF 
YYYY  < I ) =RGH 
XXXX ( I > =CNS 

IF<  I * EU . I END)  RMAXORSH 

IF<  I *E(S*  [END)  RMAX1-RSH1 

IF ( I . EQ . I END ) RHAX»2.0#RMAX t-RHAXC 

IF  ( DABS  < XNS-XNSO  > * (3  T . TF ACT  ) TFACT-  DABS  ( XNS-XNSO  ) 

IF  ( S<GE  *0.0001 ) GO  TO  641 
DO  640  N-l i IE 

RCSF ( N ) -CNS/ ( 1 *+CK*ONS#XN(N) ) 

IJ1  (N>-U2(N) 

V1(N)«U2<N) 

PI (N)=P2(N) 

Ti <N)=T2(N> 

R.l  < N ) -R2  ( N ) 

UC<N)*U2<N> 

00 ( N > =V2 ( N ) 

PC  ( N ) -P2  ( N ) 

TO ( N ) =T2  < N ) 

RC(N)=F'C(N)/TC(N) 

UIN ( M )=U2N ( N) 

U1N(N>=V2N(N> 

PLN(N)=P2N<N) 
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TIN ( N ) - T2N ( N ) 

UCN ( N)=U2N ( N ) 

UCN(N)-=U2N<N> 

PCN<N.)  =F'2N  ( N ) 

TCN ( N ) =T2N(N) 

U LNNUN>=U2NN<N) 

T L NN  < N ) » T2NN < N ) 

640  CONTINUE 

641  CONTINUE 

SOLVE  N MOMENTUM  EQUATION  ( STAGNATION  REGION  ONLY  ) 

c*#****#**##*###**#**#*******  **######***#**#**#*#**##■* 

c 

IF < S.GE.  0,0001)  GO  TO  SOO 
IF  ( I THIN , EG  » 0 > GO  TO  716 
tF < NTIME . NE . 1 ) HQ  TO  713 
UFO -VSP 
VUS2G=UVS2 
UUSG ( 1 ) =VUS 
DO  710  N-l » IE 
UG  < N ) -U2  < N > 

URS(N) -US(N) 

710  UO < N) =VC<  N ) 

GO  TO  714 

713  UUS2G- ( UVSG < 1 ) +UUSG ( 2) )/2,0 
UPG  =< USPP( 1 >+USPP(2)  )/2.0 
DO  712  N— 1 » IE 

VG<N)=(UCI (N, 1 )+UCI IN,  2) >/2.0 
UG3<N)=(UCI <N»2)-UCI <N. 1) >/DS 
712  UO  < N ) =UC I ( Nf 1 ) 

714  CONTINUE 

DO  720  N-2,IM 

UON  < N ) = ( DN  < N - 1 ) *U0 < N f 1 ) /DN ( N ) -DN ( N ) *U0  < N- 1 > /UN < M-l ) ) / < DN ( N ) + 

1 DN(N-i) )+(DNCN)-DN<N-l ) )*V0<N>/(DN(N)*0NCN-1 ) ) 

UCN ( N ) « < DN ( N- L ) *VG ( N f 1 > /DN « N ) -DN < N ) *UG ( N~  L ) /DN ( N-l ) ) / ( DN  ( N ) + 

1 DN(N-l) ) F < DN  < N ) ~DN  < N-l ) ) *UG < N > / ( DN < N ) *DN ( N- 1 ) ) 

720  CONTINUE 

UON  (IE)  =U0<  IE  )*<DN(IM-i>  F2 . *DN ( IM ) ) / < DN < I M ) * ( DN < I M ) f DN  ( L'M-D  ) > 

1 -U0<  IE-1  )*<DN<  I M-l  )+DN(IM)  ) / ( DN  < I M ) #DN  < I M -1 ) ) 

2 +U0< IE-2)*DN( IM)/(DN< IM-i >*<DN( IM> FDN< IH~1 ) ) > 

UON ( IE)=UG< IE)*CON< IM-t )+2.*DN( IM) ) / ( DN C IM ) * < DN ( IM ) FDN< IM-l > ) ) 

1 -UGI IE-1 >*(DN< IM-l )+DN( IM) >/<DN( IM)*DN< IM-l ) ) 

2 t UG  C IE-2 ) *DN ( IM)/ ( DN ( IM-l ) #1 DN  < IM  M DN( IM-1 ) ) ) 

UQN<  1>-=~U0(1  >#'DN<2)+2,*DN<1  ) )/<ONCl  )*(DN(1  ) FDN<2>  ) > 

1 +U0  ( 2 > * ( DN  < 2 ) FDN  ( 1 ) ) / < DN  < 1 ) *DN  <,  2 ) ) 

2 -UO  C3  >*DN( 1 >/<DN(2) *(DN<1 ) FDN ( 2) ) ) 

UDN ( 1 )=-UG( t)*(DN(2) F2.*DN<1> ) / < DN <1 ) * < DN ( 1 > FDN ( 2 ) ) ) 

1 -*UG<2  WDN<2)  TDN<1))/(DN(1)*DN(2)  ) 

2 -UG(3)*DN< 1)/CDN<2)*<DN< 1)+DN<2) ) ) 

716  CONTINUE 
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c 

C ******  **********************  ************************* 

c 

PJ  IN  t IE  ) "NI\lS2*UUS2**2  *CK2  *XN!i/  ( PFNS2*  < 1 , +CK2*XNS>  > 

P2U  IE> - t .0 

pan<  le> --I 'i'iN ( .n: 5 

PA<  I.LZ)  -t  , 

PON ( [0-0,0 
l-’O « F E ) - L * 0 

CAL  L SHOALS  ( 1 . 000  t 0 , ODO  > 1 ♦ 000 1 0 . 000 1 1T  CO  > 0030  < IJUSO » FT  SO  » t ) 

ir<  i ri 1 1 n ► t-ro . o > go  ro  ?o o 

F'C5.W(Li;  )- -RKS2*WS2G*0082fi*<  <R2<  IF  >*VG(IE‘  )™R2(IE  >*U2(IE  )*UUG2 

1 #XNSPM*XN ( IE  i / ( U0R2G* ( 1 . FCK2*XNS*XN (IE  ))>>*VnN<IE  ) 

2 HJUS2*XN£J#R2<  EE  )*IJ2<IE  >#(OGS(IE  >+UPG*VG< EE  )/00!3?G> 

3 / < 0US2G*  < L ♦ FCK2*XNS*XN<  IE  ))>)/FrS2 
F' .33  (IE)  =0,0 

I' AN ( [E>--PAN<IE>fP3.3N(  IE) 

PON  ( IE  ) -WOO  ( t ) *P0  ( I E ) * 0 0 ( IE  )*VONlIE  ) / < PPSO*  T2  ( IE  )) 

750  CONTINUE 
KON  = IN 

00  GLO  N ■ - 1 s>  1 M 

F''2  IN  ( NON  > -RRC2*UUC2**2  *CK2*XNS*R2  ( KON ) *02  ( NON  ) **2  /<PF:'S2*U  »+CK2* 
1 XNS*XN ( NON ) ) ) 

P21(K0N)--K2.l  ( KON  F t >-DN< KON )* < F''2  LN (KON+t ) FF'2iN(K0N  ) ) /2* 
PANtK0N>=P21N(K0N) 

PA ( NON ) -=P2  L (KON) 

PON (KON) =0.0 

[T(  ITH  IN . EG , 0 ) CO  TO  305 

P3.SN  < KON ) = ~RRS2*00!:J20*V0S2G*  ( ( R2  ( KON  ) *013  ( KON  ) -K2  ( KON  ) *U2  < KON  > *!1UB2 

1 *XNSPM*XN(  KON ) / ( 00G2G*  ( 1 . F CK2*XNS*XN  < KON ) ) ) ) *OGN  < KON ) 

2 +UUS2*XNS#R2  ( KON  ) *U2  ( KON ) * < OGS  < KON ) f OPG*OG  ( KON  ) /V0S2G ) 

3 / ( 00S2G*  ( t ,+CK2*XNS#XNCK0N)  ))  >/PPS2 

P3.3  < KON ) -F'33  < KON  f 1 ) -DN ( KON ) * ( P3  5N ( KON  H ) +P33N ( KON ) > /2 , 

PAN ( KON ) =PAN < KON ) < P33N  < KON ) 

F'A  ( KON ) =PA  < KON  ) f P33  ( KON  ) 

PON < KON > --OOCG ( 1 )*l  2 ( KON  > *00  < KON ) *UON  < KON ) / < PPSO*  T2 t KON > ) 

005  CONTINUE 

|-‘0(  KON)  =P0<  KON+  l > —UN  C KON)  * ( PON<  KON-f  1 )+F’ON  (KON ) )/2. 

0 I 0 KON- KON~ l 
000  CONTINUE 
C 

C****************************  ************************* 

c 

IF (NT TT .GT.200)  GO  TO  3000 
IF ( 1 F ACT , G T * Xf'ACT ) GO  TO  2000 
K’HFAO  "RRS*OOh*CN5/  ( EPS*EPS*V  l SCO  ) 

CEUH--?  . 0 * I J U S * R R ! i * 0 0 M * 0 E S C ( 1 ) * ( IJCN  < \ ) - CK*CNS*UC  ( t ) >/k'EKAC 
HEAr~-rUi*RRG*OOM*';CON(.IKOIGC(I ) ITCN(  t ) /CISCO  fUIJS*UUS*VI  GC  ( L )* 

L UC  < l > *UCN  ( 1 ) / T TS ) /REFAC 

STAN  =HEAT / (0,5 F l .0/ ( ( GAM-i , 0 ) *KftAO*RMAC ) - TU ) 
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IFCI.LE.2)  VVSG(I)-VVS 
IF  ( I . I.E  .2  > UGPP(T>=VSP 
CRN  I - l .orio 
Dl)  840  N = 1 1 1E 

XM  < N)-=BSGR  r ( OJUS*UIJS*UC(N)*UC(N)  FWM*OVM*VC  ( N ) *VC  ( N ) >/(  (GAM-1.0) 
1 *TTS*rC(N> ) ) 

P02P01 -1 . 0 

IF(XM(N)  *l.E*l  .0)  GO  TO  799 

P02P01  = < < GAM  FI  . 0 ) #XM  < N ) *Xh  ( N ) / ( 2 . 0 + < GAM-1 . 0 ) *XM  ( N ) *XM  ( N ) ) > ** 

1 (GAM/ (GAM-1 .0) >/(2.0*GAM*XM(N)*XM(N)/(GAM+l ,0)-<GAM-l .0)/ 

2 (GAM+1 . 0) )#*( 1 .0/<GAM-l .0) > 

799  PITO(N>-P02PQ1*PC(N)*PPS*( 1 . 0+( GAM-1 . 0 ) *XM ( N ) *XM ( N ) /2 . 0 ) ** 

1 ( GAM/ (GAM-1.0) )/P0 IP 

IFCI.EQ.l)  VC(N)=VC(N)*VVS 
IF(I.EG.i)  V2(N)=VC<N) 

IF(I.EG.l)  VI (N  >=V2 ( N ) 

IF(  I ♦ L.E  * 2 > VCI(N»I)=VC(N)/WS 
US ( N > = < U2 ( N ) -U 1 ( N ) ) /OS 
VS(N)”(V2<N)-V1(N) > /DS 
PS<N)“(P2<N>-P1 <N>)/D3 
TG ( N ) = < T2 ( N ) -T1 ( N ) > /DS 
UC ( N )— ( U2 ( N >+Ul <N) )/2»0 
UC(N) -(V2CN>+V1(N>  >/2.0 
PC ( N ) = (P2 (N ) +P1 ( N) ) /2 .0 
TC ( N ) — ( T2 ( N) +T1 ( N) ) /2  * 0 
UCN  ( N ) = ( U2N  ( N ) -t-U  1 N < N ) )/2.0 
VCN(N>=(V2N(N)+V1N<N) )/2.0 
PCN(N>-(P2N(N)  FP1N(N>  >/2.0 
TCN( N ) = ( T2N  ( N >+TlN  CN  > )/2.0 
RC  < N ) =PC(N)/TC( N) 

U.L  ( N>=U2  (N) 

VI (N)=V2(N> 

T1 (N)-T2(N) 

R1 (N)-R2(N) 

PI (N)=P2(N> 

T1N(N)=T2N(N) 

T1NN( N) -T2NN<  N ) 

IJ1 N ( N ) -U2N  ( N ) 

U INN ( N ) =U2NN ( N ) 

V1N(N)=V2N(N> 

P1N<N)“ P2N(N  > 

C01=C02 
S40  CONTINUE 

PUALL=PPS*PC ( 1 > 

IF<S. LE. 0*0001 ) GO  TO  841 

C0P2=4.O#kS#SIF*PUALL 

CDF2-2.0.*RS*CSF*CFCH 

CDPP-=CDPD  + <C0Pl+C0P2)#DS/2*0 

CDFD=CDFD+(CDFl+CDF2>*DS/2.0 

OBP=CDPD/(RS#RS> 
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CUF-CDI-  D/  < KS*Ri> ) 

(ML  IF iS.OE ►0*0001)  UU  TO  8 4? 

F'UALO-"PWALL 
CDI:  ^0*0 

onp  *2.o*puall» 

04 2 CD!  0T--CDF  fCDR 
CUM.  -C0P2 
rriFi -!.;riF2 
PWRaT  sF’WALL/F'WALO 
C 

0* *************** ************  ************************* 

C 

XNSMXNS 

CNS-XNS 

UUSI~'UUS2 

vvsmvvs? 

rfCMTrS? 

PPGt-“PPS2 
R kSi--RRS2 
RS1--RS2 
CSF L=CSF2 
OKI ®CK2 

ATi  ( I ) ~ -2  » 0*CK*DT AN  ( ALPC-PH  TO ) -CNS*CKP/  (1.0  fCK*CNS ) 

AB(I)=-2.0/Dr 

AC  a >--YNSPP< I >+<2.*CK*DTANl ALPC-PHIC) fCNS*CKP/<l ,04CK*CNS) > 

1 *YNSP < I ) K2  »0* ( 2 * Q*RSH -YNGH ( I ) V DT 

IF(NTt«E»NI£.NTIMEl>  GO  TO  986 
TF<  < I - l )/  [PRINT*IPRINT  ► NE.  ( I-  1 ) ) DO  TO  906 
XF"(J,Gr*t  l>  l PRINT -I PR  [Nt 

UR.t  TE  C A r 908  > S > XB  t RS , CNSJ,f  XH!  I » R5H  r EPS  » N T T T > N TIME 

908  FORMAT (//BX  F!i'i  10X*  'X'  * 1 OX » 'R'  >9X,  'NSH'  riiX?  ' XSH ' »0X»  ‘'ROM  ' , 

1 OX  j 'EFS'  vAX?  'NO  ITER'yAXf  ' N T1  ME  ' /.TX  > 7 ( F9  ♦ 5 » 2X ) *16*112) 

WRITE  ( A / 909  ) CFDH  * HEAT  * STAN  * CDF  * COP  i CD  TO  T , F'WALL  * F'WRA  J’ 

909  FORMAT  < /AX* 'CF' *9X» 'HEAT' * 7X , 'STAN' * 7X» 'CDF' -OX* 'COP' *8X» 

1 ' CD  T OT  ' * 5X  * ' PUALL  ' * AX  * ' PW/PO ' /3X  * 8 ( Fr9  ♦ 5 * 2X ) ) 

WRITE <6* 910)  UUS* VVS.PPSj ITS*RRS 

910  FORMA  f < /AX  * 'UUS' »8X* 'VUS' ,8X> 'PPS' *8X* 'TTS'»8X* ' RRS ' / 

1 3X*S<F9.5*2X) ) 

WRI TEC  6 *911 ) 

UiR  1 TE<  6*9  L2)  ( XN  ( N ) * UC  < N ) » VC  ( N > * PC  < N ) * TC  C N ) * RC  ( N ) » XM  ( N ) * 

1 PITO(N) *N=1 * TE) 

911  FORMAT  ( /8X  * ' N/NSH ' * 7X  * ' IJ/USH ' * 7X  * ' V ' * 7X  * ' P/PSH  ' * 7X  ? ' T/TSH ' * 

1 7X  * ' R/RSFI ' ? 7X  * 'MACH'  *8X/  'PTTU'  ) 

9 J.  2 H ORMAT  (A  X * OK  L 2 . A ) 

986  CONTINUE 
l?S"RS2 
S-F3FDS2 

CALL  GEOM  ( S . DS2  * RS  * CK  , CSF  * S TF'  > XIO 
US2“RS 

PHI  C-DARCOSCCSl-  ) 
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RSH=RS 1CNSACSF 
5000  CONTINUE 

r;  ***,*$$*'4>$$*‘t******  N+l  TIME  SWEEP  *************** 

CALL  MANISH  <BS* IEND  > RMAX » CONU ) 

IF(NTTME.NE.NTIMEt)  GO  TO  987 
WRITE (6*914) 

UR I TE ( 6 » 915 ) < YNSH  ( I > > 1 = 1 » IENLU  ) 

UR [ TE ( 6 y 91 6 > ( XNSH< I > » 1 = 1 > IENLU  ) * 

915  FORMA  TC/SXi- 'SHOCK  RADUIS ' /30 ( 3X > 10F10 . 4/ ) ) 

91-1  FORMAT  <////////» 3X * 'SHOCK  SHAPE  ' ///AX  FI  NAL  SWEEP') 

916  FORMAT T/SX* 'SHOCK  TH ICKNESS ' > 730 < 3X » 1 0F10 , 4/ ) ) 

WR [IE (A >913)  CQNUtNTXME 

913  FORMAT ( /3X> 'FINAL  SWEEP  HAS  CONVERGED  TO  : '»F12.4,2X> 

1 '(CHANGE  IN  RSH  / RSH  > f'»4Xr'N0  TIME  CYCLES  =',I8/) 

WR  I TE  ( A f 91 7 ) 

WRITE (A >915)  <YYYY(I) t 1=1 rlENIH) 

WR I TE (6*916)  (XXXX(I)  > I =1  > IEND1 ) 

917  FORMAT <//4X> 'STAR  SWEEP') 

987  CONTINUE 

IF < CONY. GT .0.001 OOOODO)  GO  TO  20 
IFINTIME.EG.NTIMEI ) GO  TO  3000 
NTIME1=NTIME+1 
CO  TO  20 
3000  CONTINUE 
STOP 
END 

SUBROUTINE  SHVALS  <SP>  CP»  SFB.  CPB»  TTSH»  ORSH r URSH  > F'PSH?  ID) 
MOD  729. It  778 >5 

IMPLICIT  REAL*8  <A-H>  D-Z) 

COMMON  /INSH/  CONO  t GAM  » S > UPSH  » XNS  > 


COMMAND?  DEL  72?.1>778.5 
COMMAND?  INSERT  728.1 >777.1 

728.1  ? t 

777.1  ? i 

COMMAND?  LIST  704/779  UNNUMBERED  MARKERS 
RSH=RS+CNS*CSF 

COMMAND'5*  LIST  704/778  UNNUMBERED  MBRKER=* 
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RSI-I -RS+CNS*nSF 
SOOO  CONTINUE 

it'Ti ds N+i  rim:  swllf'  ^ ^ *c. -*11  <»> ,-{<■  'Kj^ * 

CALL  MAN  I SI-1  ( DG  r [END » RMAX  r CONV  > 

IF(NTIME*NE.NTIMEi)  GO  TO  907 
WRITE < 6 » 9 L 4 ) 

WRITE,  I A ? 9 Li:< ) <YNSIR  D »I-i  7 I END  L ) 

UR I TE< 6*916)  < XNSH ( J 1 * T~1 » IEND L > 

9 t5  FORMAT  < /15X  •>  ' SHUCK  RATUJ I S ' /30  ( 3X  • LOLL  0.4/)) 

9U  FORMA  T < //////// * 3X  ? ' SHOCK  SHAPE ' ///AX  FINAL  SWEEP') 

91 A FORMAT (/SX 7 'SHOCK  TH1CKMEi30  ' *730(3X7  10F10*4/  > ) 

WRI TE ( A » 913 ) CONV7NTIME 

9L3  FORMAT < /OX 7 'FINAL  SWEEP  HAS  CONVERGED  TO  % ' »f L2*4 »2Xr 
1 ' { CHANGE  IN  RSH  / RSH ) *'»4X7'NO  TIME  CYCLES  ='7.[0/> 

WR I TE  ( A * 9 1 7 > 

WRirEiA,9IS>  < YYYY < I > r T-J.  * TEND1 ) 

WRITE  (A?  916)  (XXXX(I)  • T--1  * IEND1 ) 

917  FORMA  T < //4X» 'STAR  SWEEP') 

9B7  CONTINUE 

IF (CONV.GT.O *0010000110)  GO  TO  20 
IF(NTIME.EO.NTIMEi)  GO  TO  3000 
NTIMCl-NTIME+l 
GO  TO  20 
3000  CONTINUE 
STOP 
END 
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LO 


SUBROUTINE  SHUAI.S  (SP. 

CP } SPB» 

CPB.  TTSH. 

URSH. 

URSH.  PPSH 

IMPLICIT  REALMS  ( A— H ? 

0-Z) 

COMMON  /INSH/  CONO  * 

GAM  > 

S 

9 

Uf'SH 

. XNS 

EPS  * 

RMAC  > 

TPSH 

9 

UISCU 

COMMON/ OUTS H/  PPS  » 

RRS  . 

t rs 

r 

UUSl 

. UUS 

PPS1  . 

RRSI  . 

TTSl 

9 

UUS2 

. UUS  l 

PSP  . 

RRS  2 > 

TSP 

f 

USP 

, UUS2 

PPS  2 . 

RSP  » 

TTS2 

f 

UUS 

. USP 

ID  ) 


COMMON/MAT NN/CRN1'  > DS  » DN  ( t 1 1 ) » XN  < 1 1 1 >'*  IM»  IE 
GAMP  = GAM  + 1.0  DO 
GAMM  - GAM  - l.ODO 
RMACG  - RMAC  * RMAC 
E.PSQ  ■»  CPS  * EPS 
SPG  = SP  * SP 
FOGG  - 4.0P0/<GAMP*GAMP) 

DEN  » RMACG  * RMACG  * SPG 
URSH  = SP  * CP  / <SP  + EPSG  * UISCO 
TTSH  = < <URSH-CP)##2  + FOGG*GAM*SPQ  + ( 2 » ODO/GAMM-FOOQ#GAMM  > /RMACG 
! - FOGG/DEN ) *0 . 5D04SP/ < SP+EPS0*CONO#TPSH/XNS > 

PPSH  =<2.0D0*SPG  - GAMM/  ( GAMtfRMACG)  > / GAMP 
RRSH  = GAM  * PPSH  / (GAMM  * TTSH) 

URSH  a ~SP  / RRSH 
GO  TO  < 20  » 5 ) . ID 
CONTINUE 
T TS2  = TTSH 
PPS2  = PPSH 


RRS2  = RRSH 

IJUS2  = URSH  * SPB  + URSH  * CPB 

UUS2  = -URSH  * CPB  + URSH  # SPB 

IF  (S  .GE.  .0001)  GO  TO  10 

UUS1  = -UUS2 

UUS1  = UUS2 

TTS1  = TTS2 

PPS1  = PPS2 

RRSi  = RRS2 

CONTINUE 

ULJS=  < UUS2+UUS1 ) /2  . OBO 
UUS=<UVS2+UUS1  >/2.0D0 
TTS= ( T TS2+ TTSl )/2.0D0 
PPSa(PPS2+PPSl)/2.OD0 
RRS= (RRS2+RRS1 )/2.0D0 


USP  = 

< UUS2 

- IJUS1) 

/ 

DS 

USP  a 

(UUS2 

- UUS1) 

/ 

DS 

TSP  a 

< T TS2 

- TTSl) 

/ 

DS 

PSP  a 

( PPS2 

- PPS1 ) 

/ 

DS 

RSP  a 

(RRS2 

- RRSI) 

/ 

DS 

20  CONTINUE 
RETURN 


END 
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rsufckou r .r ne  oeum ( s < ns  > fts » ck - oof  . oi r . xro 

0 E OR  A HYPIIkDUl  nin  ACYhlTUTiC  TO  A CONE  OT  TOTAL  tNTEREOR  ANGLF 

C OF  45  OEGRL'F  0 

IMPLICIT  REAL  nj  ( A -H  * (J-Z) 

ANG--22 . 5D0/57  . 29'J77VSD0 
ci-n  r an  *;  anu>#*2 
02--1  .odo  h:i 

DERR  • DO}Hri!5ilRT<Cl)Klii5*HSH  J .DO)  ZOSOk T ( C2*RG*RS+ 1 . F«> ; 

RL  Zi"=R!i ! ULRR/O.ODO 

200 1 CONTINUE 

Dir  ft  ft  ■=  DS*tiSQRT(C  l*h'EXr*REXP+l  .DO)  / D!3CH<  T < G2*REXP*REXP+1 . DO  > 
DEL  T-ftS  fHh'ftR/2  . OLiO  -ft  EX  ft' 

IF (DEL T.LE.O.OOOOt)  00  TO  2002 
RFXR-RS  Hii:Rft/2.0D0 

00  TfO  700.1 

2002  kO'-KO-fOERR 

SUIT  ~ DSliR T' ( 1 . ODO  + C2  * RS  * HS  ) 

XB—  ( - 1 . ODO+DSUR  r ( L . 000-f  C t *RC*R3  > ) 7C  l 
CK-I  .000/ (SMB  I *S(.lftr*SURT) 

CKr-kG/SORT 

<>rr  - DSUkm.ooo  - csr  * csf> 

RETURN 

1 .Nil 
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subroutine  w iwn  < tfact*x  ) 

LhPL  IfJLT  Ri:AL#G<A-H*0-Z> 

LGMhON/DEQSG/  U2  < 1 1 1 ) * V2  < LI  1 ) *P2dli  > *T2d  L L ) *U2N<  1 1 L ) * U2N<  111)* 

L P2N<  ill  ) *T2N(U1)  *U2NN<11  l ) *T2NN<  111  ) 

COMMON /BENS/A  11  d 1 1 ) * A13  d l L ) * FU  1 d 1 1 ) * B L2 d 1 1 ) » I<  13  <1 1 1 > * DI  4 (1  L L > * 

1 C l l ( 1 l l ) *C  L3d  1 1 > * [I  L < It  l ) *A21  ( 111 ) * A 2 3 ( 1 1 l ) *A24(  1 1 L ) * DJI  ( 1 l 1 > * 

2 B22  <1  1 1 ) * 023  < 1 i l ) * B24  Cl .1 1 > * C21  ( 1 1 L ) * 023  < 1 1 L ) * 024  < L 1 I ) * 02  1 1 J L ) « 

,3  031(1  LI)  *632(111)  *033(11  l ) *034(111  ) *031  ( 1 11 ) * 032  < 111  > * 033 ( LI  I ) * 

4 C3  4 ( 1 LI  > * DJ  < 1 11  ) * A 4 i C i l L ) * A 42  ( L 1 1 ) •*  A 4 3 < 1 L 1 ) *'.VH  (L  l L ) * 04  i ( L 1 L ) » 

5 B42C1 ll ) *B43( 1 LI ) *U44 (1 L 1 > * 04 < L 1 L ) 

OOhhUN/MAINN/CRN [ * D5  * DN (111) * XN < 1 1 1 > , IM * IE 
COMMON/ BU UN NO/  DU ( 1 1 1 ) * DU ' 1 1 1 ) > OP ( 1 1 L ) * 0 T < 1 1 1 ) » 

1 HU ( 1 3 1 ) * EU< 1 U ) > EP  < 1 1 1 > * E F ( 1 1 1 ) » 

FUdll  ) *FU<  l LL ) * FPC111 ) *FT<  111 ) > 

GUdll  > *GU<111 ) *GP(11  L ) *GTdl  1 > * 

HUd  Li)  *HVd  PL ) *HP(1  LI)  * HI  (lid  ) 

N-l  E 
Nl-N- L 
OUt N) -0.000 
DKNJ-O.ODO 
FLJ  ( N > ~0 . 000 
Oil  ( N )-0 . 000 
HU<N)=U?(N) 

OU < N ) - -A4 1 < N > /B42 < N ) 

EU(N)  = -A42 ( N ) /042  ( N ) 

FU ( N ) --A43 ( M ) /B42 ( N ) 

GU(N)--~A44(N)/B42(N) 

HV<N> -<04<N)-B41  (N  > *IJ2  CN ) -B43 < N > #P2 < N ) -B44 < N ) * F2<  N)  )/042(N) 

DP(N)=0.00 

EP ( N ) =0 . 00 

FP ( N ) -0 . 00 

GP ( N ) -0 . 00 

HP ( N ) -P2 (N ) 

OT(N) =0.000 
ET(N  >=0.0I>0 
F f ( N ) -0 « 000  . 

GT(N)=O.ODO 
HF(N) -r2(N) 
n ACT -0 . 000 
P20=-P2  < 1 ) 

00  30  J 1 = 1 * N 1 
J=N~J  L 

B31B=B31( J)+C31< J)*UU< J+l)+C32< J)*DV< J+l >+033( J)*DP( J+  L ) FC34( J)* 

1 DT ( J+l ) 

B32B -032  ( J ) fC31(  J5*EU<  J+l)  1032  ( J ) #EV<  JU  ) f 033  ( J ) *EP  (J+l  ) 1034  C J>* 
t ET(JH) 

B330--B33  ( J ) 1031  ( J ) -VFU(  J1  1 ) +C32  < J ) *FU  ( J FI  ) +C33  < J ) *FP  ( J H > tC34<  J)  * 

.1  F T ( J 1 L ) 

B34B=B34(  J) 1031  ( J > *GU  < JU  >+C32<  J+i  ) 1033 ( J > *GP ( J 1 1 ) 103  4 ( J)  ■!< 

1 GT(Jtl) 
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D3B  -=D3  < J ) - C3  L (.  ) -CS2  ( J ) *1 IV  < J+  L >-C3  J<  J)  4!  J|-*<  J4  1 > --C.-5-H  J ) * 

1 HKJ+I) 

.1  j" t j » lu  . n nn  to  20 

«j  iB-sui j) ft: i n j>  kdu< jh. :■  tc 13' jn-irt  ju> 

B1.3B=B  t 2<  J)  MOL  L ( J)  XEU  UH  M0L3l  n fit  Pt,J+  I J 
D I SB  =0  I ,?(  J)  ICI  I ( J)  tFlJ  Utl  ) H,  1 0 < J ; -SEP  (JU) 

0 | 4B  ~h  L 4 ( J ) f C 1 1 ' J j *GU  J K ) f C I 3 ( ..P  KHP  < J4  1 > 

(I  t 0 =D1(J)  -L‘LH  J)*I!U(JH  > -01 3( Jft > 

B21B-R21  ( J)  fC2JU  J)#DUt J+  L ) f023  < J > *RP  < J f 1 > +034 < J ) *1)  T < J M > 

B22B— D22  ( J 1 +C21  <J1#EU<  JM  >+C23<  J>*EP(  JH  >+C24<  J>*r.r<  .J  + t > 
B23B-R23<  J)  Kt21<  J5+FUCJ+1)  +C23 ( J > *FP ( J H ) H:24<  J>*rt  ( J+i) 

B24B-B24  ( J)tC21(J)  4GI.J  ( J f 1 ) f C23  < J 1 *0P  (JH>  +02  4 (J)*(iHJ+l) 

B2B  -B2  ( J ) -C2 1 ( J ) 4HU  ( J+l ) -C23  < J ) *HP  < J+L ) -024  ( J ) *M‘l  UH) 

XK1 =0330*044 < J)  -H34B*B43< J> 

XK2 =B32B*B44< J) -B34B*B42< J) 

XK3  -B32B* B4 3 < J > -B33B*H  42  ( J ) 

XK4-B23BSB44 ( J ) -B24B*B43 < J) 

XKG-B22B4044  ( J)  ~B24B#B42’<  J) 

XK6~R22B*B  43  < J ) -B23B  + B42  ( J > 

XK  7-B2.3B#  B34B-B2  IB  *B33B 
X K 8 - B 2 2 B * B 3 4 B - B 2 4 B * B 3 2 £< 

X K V =-  B 2 2 D 4t  B 3 3 D ■ - B 2 . 5 U * 6 3 2 w 

XK 1 0 -02  LB *044 ( J) -034 R# B4  L ( J ) 

XMl*R:HB*B43<J>-B33B*84i<J> 

XK12-B2  LB#R44<  J>-B24B*B4  I.  < J) 

XK.t  3--B21B&B43  < J)  -B23BKB 41  ( J ) 

XK14-»B21B#B34B”B24H*B31B 
XK1 5-B2J  BW33B  -B23B*B31  B 
XKL6=B3IB#B42( J)-B32B*D41( J) 

XK  J 7 =B2 LB#B42< ,J>-B22B*B4i  < J) 

XK  L0- B210*B32B  -B22B*B31B 

Bt JC-B22B*XK1  B23BtXK2  + B24L«#XK  J 

B21C-B  L2B*XM  -HI  3B*XK2HU4B*XK3 

03 l C =B I 20 XX K 4 B i 3B*XK5+B 1 4H*X1\6 

B41C-BJ  2B*XK7  ~H13B*Xk8 t-BJ.  4R*XK9 

BE  T -0  I.  I B*B1  LO-D21B*B2JC+B3J.B*B3JC-B4l<J)*04lC 

M1.C-B11C 

r r2 1 C ==  0 1 2 B X X K J - B 1. 3 0 X X K 2 + B .1 4 B * X K 3 
H3LC=Hl2B*XK4-B13B*XK5+B14B4tXK6 
04  tC-B120N<XK7  - B13BXXK8  HB1 4B*XK9 
D1 2C-B21  BXXK.l -B23BXXK1 0 LB24BXXK i 1 
D22C~Bi iBXXK  L -01 3BXXK 10  f B1 4BXXK1 1 
D32C-B J iB*XK4 -B J 3BXXK12 fB14PXXK I 3 
D42C-0I  1 B*XI\7-R13R*XKi41-B1  4B*XKJ5 
D L3C-B21 B»XK2-B22B+XK10 IB24BXXK 16 
D23C=B !,  t.B*XK2-B12BXXK10  + B l.4B*XK16 
0330 -B  ! I BXXK'.-j- -01  2BXXK  L2  \ l<  I 4BXXM  7 
043C  aRJ  LB#XKS  -Bt  2D  HXK14  H3  I 4B*XKtf? 

D J.  4 C ~ 0 2 1 B X X K 3 - 0 2 2 D X X K ,1  HR23B*XK  13 
R24C-H  l tU*XK3-B12U*XKJ  l + B13FUXMi4 
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D34C=B1J  B*XK6  -0l2B$XKi3  fB13EtHXKi7 

li44C-imB*XK9-Bt2B*XN15+B13P*XKl3 

HU < J ) = ( -A. I 1 < J ) #Dl  1C  f A2  l ( J > $ 021 C IA4 I < J) *D4 1C ) /BE  T 

E U < J > “ A 4 2 < J > + D 4 .1  C /D  E T 

FU< J>=<A23( J)#D21CfA43< J>*D41C-A13< J)*DJ  .LO/nET 
GLK  J>=*<  A24<  J>#D21CfA44<  J)*D4  1C) /BET 
HU(  J>  = < Li  I B*D1  I C-D2B*B2tC+B3B*D3i  C-D4  ( J ) *04 1 C ) /BET 
0U( J)=< Ail< J)*D12C-A2t ( J)*B22C-A41< J)*D42C>/DET 
EU<  J)--A42<  J>.*B42C/0ET 

FU  < J > =-  < A23  < J i *D22C  1A43  ( J ) JKD42C-A  J 3 < J) #D12C ) /BET 
GU  < J ) -- < A24 ( J ) *B22C+A44 ( J ) *D42C ) /BET 
HU  < J ) --  ( -D 1B*012C+02B*D22C-B3B*D32C+D4  < J ) *D42C )/DET 
DP<  J)  = (-A1 1 ( J)*D13C+A21  ( J)*D23C  + A4:l  < J>*D43C)/DET 
EP  < J ) =A42 ( J ) *D43C/BET 

FP< J)«(A23( J)*D23CfA43( J)#B43C-A13< J)*Di3C>/DET 
GPt J>=<A24( J)#B23C+A44( J>  #0430 /BET 
HPC  J)-<D1B*D  L3C-D2B*D23C+D3B#D33C-D4<  J >*0430 /BET 
BT  < J ) = < A t 1 < J)*B14C“A21 ( J >#D24C-A4 l ( J ) *D44C ) /BE  T 
£T<  J)  = -A42(  J>*D44C/DET 

FT  C J )=- ( A23 < J ) $024 C tA43 ( J ) *B44C-A13 C J ) *01 4C ) /BET 
G f < J ) -- < A24 ( J ) *D24C I A 44 ( J ) *D44C ) /BET 
HT< J)- 1-010*0140  K)2B*D2  4C-D3B#D34C+D4(J>*D44C>/DEr 
GO  TO  .30 

20  l-,2(l)  = <D3B-B31B*  U2<1)~  B32B*  U2  < 1 ) -B34B#T2  ( 1 ) ) /B33B 
30  CONTINUE 

TFACT = DABS < P2<l>-P20  ) 

DO  40  J=2  s N 
U20-LJ2(  J> 

U20=U2 ( J > 

P20=P2  < J ) 

T20=T2< J) 

U2(  J)=DU<  J)*U2<  J-1)+EU(  J)*U2<  J-D+FUC  J>*P2<  J-1)+GU(  J)*T2<  J-l  ) + 
1 HU ( J ) 

U2(  J)-=DU(  J)*U2t  J-1)+EU<  J)*U2<  J-1)+FV<  J)#P2<  J-l  )+GU<  J)*T2<  J-l ) + 
1 HU ( J ) 

F2< J)=0PI J)*U2< J-l>  J-EP<  J ) #U2  ( J-l ) TFP  ( J)  *P2  ( J-l  > +GP  < J > * T2 ( J-l > f 
1 HP ( J ) 

T2(  J>=DT<  J>*U2(  J-l  > fET<  J>*V2<  J-l > +FT < J ) *P2 < J-  L > +GT  < J > *T2  < J-l ) + 
1 HT  < J ) 

IF  < DABS (U2 ( J > -U20) » GT . TFACT  ) TFACT=DABS < U2 < J)-U20) 
IF(DABS(U2<J)-U20) *GT.  TFACT  > TFAC T=DABSI U2 < J >-U2U ) 

IF  < DABS<  F2  ( J )-F'20)  * GT « TFACT  ) IFACT =DABS  ( P2  ( J > -P2G  ) 

IF <DABS<T2<J )-T20) . GT • TFACT  ) TFACT  = DABS < T2 ( J ) -T20 ) 

40  CONTINUE 

HO  50  N=2  > I M 

U2NNCN >=2. 000*<U2<N+1 ) /BN < N > +U2 < N-l )/DN(N-l > )/<DN<N> t-DN  < N-l ) ) 

1 -2  <■  ODO#U2  i N ) / ( BN  ( N ) #DN  ( N-.L  ) ) 

T2NNCN)=2.0D0>K<T2<N+l)/riN<N)  + r2<N-l>/BNCN-l ) ) / ( ON  < N)  +0N  ( N-l ) ) 

.1  -2.0D0*T2CN>/<BNttmDN<N-l  > > 

U2N<N>-< DNCN-l )*U2(N+1 )/BNt N> -DN (N )$U2 ( N-i )/DN(N-l ) )/(BN<N)+ 
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t DN(N-l ) ) FCCNCN)  -DN(N-1 ) ) *U2 ( N) / ( UN (N> #DN < N-.t  > > 

02N<  M)-( DN<N-1 )+02<N+l ) /DN ( N ) -DN ( N ) *02 < N -1 ) /DN  < N - 1 ) > / < BN  < N > F 
I DN ( N-l ) ) + (DN ( N ) -DN < N-l ) ) *02< N ) / ( DN < N) *DN< N-t ) > 

P2N  ( N ) - ( DN  < N-l ) *P2  < N+l  ) /DN  ( N ) -DN  ( N ) *F‘2  (N-l  ) /DN  (N-i>> / ( DN ( N > F 
J.  DN  < N-.L  > ) + ( DN  ( N > -DN  ( N-l  ) > *P2  ( N ) / < DN  t N ) *DN  ( N-  1 > ) 

T2N  < N > = ( DN < N-l ) * T'2  ( N + l ) /DN  ( N ) -DN  ( N > * T2 < N-l  ) /DN  ( N-l  > ) / ( DN  < N ) + 
l DN  < N-i  ) ) h (DN(N)--DN(N-l ) >*T2(N)/(BN(N>#DN  (N-l ) ) 

SO  CONTINUE 

U2N(.Ll---U2(l>*<riN<2)  + 2.0D0*IiN(i  ) ) / ( DN  < 1)  * ( DN  ( 2 ) +HN  ( 1 ) ) ) 

1 FU2(2)*(DN<2>+DN( 1 ) )/ ( DN < 2 ) *DN ( 1 ) > 

2 -U2(3)*l'iN(l)/<DN(2>*(DN(l  >+DN(2)  ) ) 

02N( l)=-02( 1 >*(DN(2) F2.0D0#BN< 1 ) > / ( DN ( 1 > * ( DN ( 2 ) +DN  ( L ) ) > 

1 FV2(2)*(DN(2)+DN( 1 > )/  ( DN  (2)  *BN  ( 1 > > 

2 -02  < 3 ) *DN  ( 1 ) / ( DN  ( 2 ) * ( DN  ( 1)  +DN  ( 2 ) ) ) 

P2N(1)=-F'2<  1)*(DN(2)+2.0D0*DN<  1 ) ) /<  DN  ( 1)  # ( DN(  2 >+BN  ( 1 ) > > 

1 +P2(2)*(DN(2)  FDN(l)  >/< DN ( 2 ) *DN ( 1 > ) 

2 -P2 ( 3 ) #DN ( 1 ) / ( DN  < 2 > * ( DN (t ) FDN ( 2 ) ) ) 

T 2N (1 ) = - T2  CL ) # < DN ( 2 ) +2 , ODO*DN (1 ) ) / ( DN ( 1 ) * C DN ( 2 > +DN ( 1 > > ) 

1 +T2(2)*(DN(2) FDN  < 1 ) ) / ( DN  < 2 ) #DN ( 1 ) ) 

2 -T2(3)*DN(  1 )/(DN(2)*(DN(l  )+DN(2)  ) ) 

IJ2N(IE)=U2(IE>*(DN(  IM-l  >+2 « ODO*ON( IM)  ) / ( DN  ( I M ) * ( DN  ( .1 M ) FDN (IM-1 ) ) ) 

1 -U2(IE-1)*(DN(1M-1)+DN(IM) ) / ( DN ( I M ) *DN < IM-1 > ) 

2 +U2( IE-2) #DN(IM>/( DN< IM-1 )*(DN(IM) FDN ( IM-1 ) ) ) 

02N( IE)=U2(IE)*(DN(IM-1 >+2.0DO*DN( IM) >/(DN( IM)*(DN( IM)+DN<  TM-1 > ) ) 

1 -02 ( IE-1 )#(DN( IM-1 )+DN(IM))/(DN(IM)*DN< IM-1) > 

2 +02( IE-2  >*DN( IM  >/ ( BN ( IM-1 ) *<DN< IM> FDNCIM-1 ) ) ) 

F-2N  ( IE  ) =F2(  IE)*(DN(IM-1  > + 2 , QDO*DN  ( Ih > > / ( DN  ( IM  ) * < DN  < IM ) FDN  ( IM- 1 ) > ) 

1 -P2(IE-l)*(DN(IM-i)+DN(IM> >/(DN(  1M)*DNUM-1  > ) 

2 +P2 ( IE-2) *DN ( IM)/(DN( IM— 1 )*(ON(IM) FDN ( IM-t ) ) ) 

T2N(  IE  )■=!£(  IE) #<DN(  IM  -1)  F2.0D0*DN(  IM)  )/(DN<  IM>*<DN(  TM)+DN(IM-1 > ) ) 

1 -T2< IE-1 )*<DN( IM-1 >+DN(IM) )/< DN ( IM) *DN ( IM-1 ) ) 

2 +T2( IE~2)*DN( IM)/(DN( IM-1 )*<DN( IM) FDN < IM-1 ) ) ) 
U2NN(i)=-U2N(l)*(DN<2>+2.D0#BN(l> >/< DN ( 1 )*< DN<2) +DN( 1 ) ) ) 

1 FU2N (2)#<  DN ( 2 ) + DN( 1 ) >/( DN( 2) *DN C 1 ) ) 

2 -U2N(3)*DN<1)/<DN(2)*(DN(1)+DN(2) ) ) 

T2NN( 1 )=- I 2N ( 1 ) *( ON (2)+2>D0*DN< 1 ) )/(DN( 1 >*< DN(2)+DN< 1 ) ) ) 

1 +T2N(2)*( DN(2)+DN( 1 ) ) /( DN< 2 ) *DN ( 1 ) ) 

2 - T2N( 3) *DN< 1 ) / ( DN( 2 ) # ( DN( 1 ) FDN (2) ) ) 

U2NN ( IE  > -U2N( IE  > *( DN ( IM-1 )+2 . DO*DN ( IM ) ) / ( DN ( IM ) * ( DN ( IM)+DN( IM-1 ) ) ) 

1 -U2NC IE-1 >*<DN( IM-1) FDN ( IM) ) / < DN ( I M ) *DN < I M- 1 ) ) 

2 +U2NC IE-2>*BN( IM)/(DN( TM-1 >*(l!N( IM) FDN( IM-1 ) ) ) 

T2NN(  IE  > ’=T2N  ( IE  )*  (DN(  IM-1  )+2*  DO*DN(  IM  ) ) / (DN(  Ih  ) * ( DN  ( IM  > FDN  ( IM-I  ) > ) 

1 ~T2N(IE-1>*CPN(  IM-.L)+DN(IM)  ) / ( DN C IM ) *DN < IM- t > > 

2 +T2N( IE-2) *DN(IM)/(DN (IM-1 >*(DN( IM) FDN(  CM-1 ) ) ) 

RETURN 

END 
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SUBROU  TINE  DEklU  CDS  r I END , I Ii ) 

CMPL1CIT  KEAL*0  <A-HjO-2) 

C0t“iM0N/L5HCI\G/XN5l  I *.  202 ) » XNSP  ( 202  ) i XNSPT-  ( 20?  > r YNSHC202)  » YNSP<?02>  y 
l YNSPP  ( 202  > » A L ( 20;?  > t A2  < 202 ) * A3  (202) 

IFNDJ.  --IEND  + 1 

xs=o  • 0 

RS-O.ODO 

DO  10  C=2fIENDl 

CALL  GEDMc  <£ y DS » KS > CK r CSF r G [ F « XB) 

IF  < ID . EQ . 1 > YNSH  ( I ) =RS  f XNSH  ( I ) *CSF 
CF  ( III  4 EH  . 2 ) XNSH  (!)-=<  YNSH  ( I ) -RS  ) /CSF 

xs=xs+ns 

10  CONTINUE 

[F<ID.E0.2>  XNSH< t > = < 4 . 0*XNSH ( 2 > -XNSH < 3 ) >/3.0 
XNSP (1 >=0.0 

YNSPF'  ( 1 ) = < 2 ♦ 0#YNSH  ( 1 ) -5 . 0*YNSH  (2 ) +4  . O&YNSH  ( 3 > - YNSH  ( A > )/<DS*l>S> 

XNSP  I-  C I.  >=<2.0*XNSH<  1 ) -5 , 0#XNSH  ( 2 >1 A . OiKXNSH  ( 3 > -XNSH  < 4 > >/(  DS#nS> 

YNSP  < 1 ) = ( 4. 0*YNSH(2>-YNSH( 3 ) -3 . 0* YNSH ( 1 ) >/<2.0*BS) 

DO  20  I =2 y I END 

XNSP (I >-=<XNSH< I M >-XNSHCI-l ) >/(2.0#DS> 

YNSP  (I  )=<YNSH<  IH>-YNSH(I~t  > > / < 2 . 0*BS > 

XNSPP  ( I ) = ( XNSH ( I+t  >-2.0*XNSHt  I)+XNSH<  I--  L ) )/<0S*DS) 

YNSPF-  ( I >--<YNSH(  1 + 1 > -2.0* YNSH < I >+YNSH<  I-I  > )/<DS*DG> 

20  CONTINUE 

XNSPt  IEND1  )=(3.0*XNSH(IEND1  )-4.0*XNSH(  I END  > +XNSH  < IEND-1  ) >/<2.0*DS> 
YNSPCIENDL ) = < 3 . 0*YN5H ( IEND1 >-4.0*YNSH< [END >+ YNSH ( IEND-1 > >/(2.0*BS) 
XNSPP < IENDl)=<2.0*XNSH( L END 1 > -5 . OYXNSH (TEND > M . 0*XNSH ( t END 1-2)- 
1 XNSH (IEND 1-3) )/(DS*DS) 

YNSPF  < IEND1 ) = < 2 . 0*YNSH  ( IEND1 ) -5 . OYYNSH  ( CENfi)+4 .0*YNSH<  I END  1-2)  - 
1 YMSHC  IEND  1-3)  )/<DS#nS> 

RETURN 

END 
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A EOC-TR  -79-25 


!iui.<kuu r i mi.,  man  i:sh  < hh ? ue  wn <■  rnax , cumu  > 
r hi 'L  Li:  i r real.**  (n  -h - fj~  n 


CO  Mil  UN/ ill  !Ck(./  XNU! ! i 202 ) > XNriP  ( ?0  2 > • XNUI  -p  ( 202)  y YMUH  ( 20, > > 


YNSPP  (202  ) y A I < 202  > > A 2 ( 202  202  ; 


ri  1 ML' Nil  L ON  E 1 202  ) . F ( 202  > 


rh  = ii£Mii 
E <.  I ' - 0.0 
1-  K t /=0 « 0 
Litl  JO  1<NTH 

A---  1 . 0/  < DS.KOS  > -A  I t t ) / ( 2 . OyDU  ) 


S--2.0/UVJ  + P!;)  f A 2 ( 3!  V 
C - 1.  ♦ 0/  ( liii'KDB  > f A U I ) / < a » OYLiO  ) 

I.  ) 

hi  ) =-c/o^ru-i ) > 

r t [ > = ■.  n-A *r  ( ] -i  > > / 1 n f a ke  u - n > 

10  conitnul 

CuMO-'O.O 
MJN= 1 M 

YNSIU  fENO  f l >=RHf-lX 
00  20  C - I.  » [ h 
YO  - i'NSH  ( klJM  > 

i'NSH  ( MJH ) ■■•IE  ( KUN ) *YNSH  < MJN+t ) fF MIN  > 
CONV  I -It  A US  < r NON  ( KUN ) -YO ) 

IF ( I 4 M E * 1 M > CONU L-CGNVt /YO 

i r t conu i .or* cunv > cono  =convi 

20  NON- KON  - I 

CAM  DLI<  LO  (Di3y  LENTt?  2) 

RI  TURN 
ENU 


y YNSP(202)  v 
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A EDC-TR-79-2S 


LU.OCK  n a r A 

IMPLICIT  REAL #8  <A-H»0-Z> 


C0MrlGN/SHt:i\G/XNSH<20?>  »XHSh,<20L‘)  >XNSPP t 202 ) rYNSM<?G:>'  t YNSf ' ( 2012 ) 


1.  YNSF'P  ( 202  ) > AIK  202  > * ABC  202  7 j AC  (202) 

HAT  A XNSH/202*0.10  72D0/ 


l£Nl.r 
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AEDC-TR-79-25 


INPUT  DAT rt 


21.7U0  0,050  430 . 000  35I.U00 

t 5 L .71  7.  J 

O .100  10. 000  1 , 400  0 , 700  O . 1 00 D -Opt 


Program  Output: 

In  the  following  pages,  parts  of  the  output  are  presented. 
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SHOCK  LAVED  PROuRAH 


1NPL1  OA1* 

PINE  TINE KtriR  TS7TC  COT  P* 

21.79  351.80  6JC.C0  li.OS  I. A O.IO 

FULL-SHCcK  LAVtK  AOi  "Cf  STEPS"  IA~N"*  51 M-.  LF  STEPS~nN~5~»  21 "S  STEP  SIZE  *0.100  “ T STEP  SIZE  - 

IMIlaL  ShtCK  SMAPc 


0^1072  C 0?it>72  0.1072  U .1072  0.1012  0.1072  0.1072  0.1C/2  0.1072  0.1072 

3.  1C 72  0.1072  - 0. 1(72 CTTC72 OnC12~  0TI0T2  0T1072  0.1072  011072  0.1072  " 

0.1072  0.11)72 


O^J  J.llOS  0.2196  0.3261  0.6251  0.5282  0.6236  0.7166  0.0021  0.8862 

0.9672  1.06S3  ~ 1.12CS T.ITU  1V7ET2  I77W  ITAOIO  1.667E  113322  1.5953 

1.4572  1.7114 


SHUCK  SEAPE 


O 

00 


FIAAl  SaEEP 


SHUCK  RALOIS 

0. 9  J . 1125 

1.417.  1.1C60 

1. elbZ  l.b908 


0.2239  0.3332  0.6397  0.5632  0.6636  0.7610  O.S35i>  0.9276 

1.1901  1.273b  1.3558  1.6Jb3  1.5153  1.5929  1.6693  1.7666 


SHUCK  THICKNESS 

J • 126 . 0.1275 

0.1293- 

0.1326 

— OrlieE — 

u.nz2 — 

— 0;l6B8  - 

0.1565 

0. 1651 

•J.l84b 

3.1457 

0.2072 

0.2191 

0.2316 

3.2661 

0.2570 

C.27UO 

0.2832 

0. JC96 

(3.3222 



' 



— . , 

EIKAL  suttP  Has  CoNVEROEO  It  » C .U009  (CHANCE  IN  RSH  / RSH1  • NC  11NE  CYCLES  « 


0.1766 

0.2966 


18 


STAR  SALLE 


SHUCK  RAbUlS 

■J.  J -3.1137 

C...03 

0.3265 

1.2763 

0.6366 

3.53E2 

0.6392 

0.7375 

1.5926 

0.B33I 

0.9260 

1.7633 

1.C166  1 . i'J65 

1 .el  73  U.O 

1.1  506 

1 .3562  ' 

1 «*i)65 

SHOCK  T.E  ICKI.-SS 
O.lOiid 

>.1836  3.1935 

0.1112 

u .20  75 

0.1162 

0.2157 

0.1227 

C.2J20 

0.1308 

2.2i63 

0. 1601 
0.2567 

0.1501 

0.2692 

0.  160b 
0.2014 

0.1720 
0 .2948 

J.3C(J.  '3.0 

10.3 


S2-6Z-UJ.-0O3V 


'*  rj  c 


5 

A 

K KSH 

ASH 

HSH 

tPS 

Nil  ITEK 

. J'Ji'O  J l 

.12054  1 

.5)168  0.31496  0.54020  1.8173) 

U. 22339 

6 

IP 

HEAT 

SIAN  COE 

CCP 

tom 

P8  ALL 

P8/P0 

.14601  0 

• 05-.J5  V 

.1120  C.  15744  1.03505  1.1924d 

0.35755 

0.36234 

005 

V8S 

PPS  IIS 

PBS 

.72558  -0 

.03027  0 

.40534  0.24158  5.67244 

N/Nih 

o/oSh 

V 

P/PSH 

1/1SH 

K/KSH 

MALII 

0.0 

J.) 

0.0 

0.86*110 

0.133039 

8.561707 

0.0 

u.02  K"30 

0.  2 6C3  74 

0.900011 

0.8821b* — 

0*22844) 

3.8b 1623 

C.  395830 

0 . 090)  Jl/ 

•3.14  1340 

-0.3)6037 

0.ee*4co 

0.317555 

2.7752)9 

0.585749 

0 • Cu  JL JO 

0.1536.19 

-0.003150 

0.86*817 

0.391015 

2.257756 

0.7*8196 

0.06CCJO 

3.240453 

-3.0C0J23-- 

”0.083399 

— C.45J667 

“1*941240 

0.839)30 

' • l J 0601 

-3.283393 

-J.OC0551 

0.884136 

C. 598 763 

1.737747 

0.9)380  7 

>.122000 

J. J 23283 

-3.0UCU3O 

0.885020 

0.556C  Id 

1.58600b 

1.C16671 

v. 140000 

0.369O63 

-0.0011 t? 

-0*86t>t)42 v*  602442 

1*470750 

1.091513 

w.160309 

0.395696 

-C. Of)  15*7 

0.68  7195 

C. 642601 

1.380202 

1.159609 

3.10000 

). 41924*. 

-C. 0)1937 

0.888474 

C. 679644 

1.307265 

1. *2*415 

U *29Ju Jj 

U.4CC694 

-C.  002364  — 

0.885872 — 

— Or7 1)391 

1*247384 

1.2808 19 

0.22  0600 

0.491)01 

-0.002865 

0.691365 

C. 744377 

1.141491 

1.335470 

C . 24  )oO  ) 

0.515662 

-0.003317 

0. B930C6 

C. 77*873 

1.155437 

1.386860 

v . 2oi>U)  J 

•J. 547929 

-0.0u)5  19 

”0**9 47)t 

~C* 1541*5 

“1711566  6 

1.415371- 

0.2 60 GOO 

0.513123 

-0.004487 

0.896555 

0.62)264 

1.089325 

1.481306 

U .339005 

•},i  5 8030 

-0.0050  79 

0.896473 

C. 645)12 

1.06*637 

1.5*4910 

9 . 3*0.)  >0 

0.6*18)8 

3*03565) 

0*900483 

—9*665554 

171)39823 

1.566387 

0.34 Of J) 

0.6  443)6 

-3.006)30 

0.902573 

C.  684034 

1.920046 

1.605507 

U .360003 

9.0661 70 

-0.0065t5 

C. 504741 

0.9C2144 

1.002885 

1.643615 

u .360000 

0.686640 

—0.0076  56 

0*406941 

—679169*3 

-0*481990  — 

1*675637  — 

C .400003 

0. 726555 

-C. 008344 

0.9093*0 

0.532560 

0.975979 

1.714082 

0.4200)0 

■ 4 i.  #.  t \ 4 .416 

0.725)44 

-0.005045 

0.511711 

C.  545638 

C.  963919 

1.74  7051 

L . H H " J JJ 

u.  / 4)257' 

’ — -*■  0 *00  97  59— 

“OlltTlA  7 

— 

0.954316 

1.778630  — 

9 . 4o  OOO'J 

3.760)09 

-C. 010464 

0.916683 

C. 568933 

0.946105 

1.808902 

0.430000 

0.7763)7 

-C. 011*18 

0.515256 

0.578620 

0.939149 

1.837940 

C.5I3L  1)9 

0. 751472 

-9.011561 

0.9*1886 

C.  587/40 

0.9333*9 

1.665dl5 

9 • 52  JC  JO 

3. 8 06642” 

-C * 0 127  T I 

8.9*4)66 

C75557 1) 

”0  .927545 

~1 .89259 1 

3.540c0‘) 

->.6*0)76 

-0.013467 

0.5*7294 

1.0C2755 

0.524709 

1.918330 

J • 56>c ) J 

*.833603 

-0.014*28 

0.530067 

1.007027 

0.921746 

1.943092 

2 .5601  Jo 

9.8  #6352 

-). 914552 

0.53*863 

1.014454 

0.915591 

1.966932 

C . 6 JC  00 9 

0.856251 

-0.01)758 

0.535740 

1.019118 

0.918187 

1.965905 

3.62COOO 

0.865527 

-9. Olo  r *3 

“075266)6 

1.1*3155 

1.41 #48) 

“2.01*064  — 

<’.640000 

0.880207 

-0.017*53 

0.541565 

1.026304 

0.91143b 

2.033456 

c.oguouj 

O.a9C)20 

-0.3 1 60  59 

0.944537 

1 .026699 

0.916007 

2.054136 

0 .6  JCOJ'J 

0.655392 

—0.0186c) 

”075475<C 

1.0)6  0 

”07915163  

2.074146 

J. 70 3000 

O.9Cd450 

—9 .3 19>64 

0.550575 

1.022256 

C.92C872 

2.09)534 

..72  JC'JO 

0.517621 

-0.029342 

0.953644 

1.033083 

0.923108 

2.11*344 

) . 74  )L  3 * 

9.925633 

— ■=). 921954 

0755T744 

1.03)72 

”0.925846 

2.13C619 

6.  76  06  )0 

0.933303 

-3.02184* 

0.555875 

1 .03) iol 

0.9*9067 

2.146399 

0.7M990J 

LI . 9 4*')6q 

-0.022563 

9.56)038 

1.C324  72 

0.932750 

*.165727 

v . 89  li'OO 

0.9  4 7442 

-3.0*3? 1 7 

“0.567  233“— 

““ 1.021232 

”0.936878 

2.18*639 

J . 82 JC  Jj 

3.95)957 

-0.0*4945  , 

0.969458 

1 .0*9  7t4 

0.941437 

2.195173 

* . b t )liji> 

C.66L  iJ3 

— w .0247 65 

L. 97*716 

1.C27793 

U. 946412 

2.215364 

i.eojro) 

0.965494 

* -9.925578 

537576005  — 

l. 0*5441 

”0795179  2 — 

2.231248 

J.6610JJ 

3.5)1363 

-). 0*6163 

0.979)29 

1.C22729 

0.957565 

2. *46656 

0.9)30  10 

0.9 76860 

- J .0*6861 

o.suioet 

1 . C 156  78 

0.463722 

2.26*221 

■..4290) 

). 961436 

-9.02  7571  

“0.506077 

[797  6208 

07970254 

2.277371 

C .94  JC)J 

I.9C6  766 

-).C*e*S4 

0.569503 

1.312638 

0.977154 

2.29*337 

./  . 5)  o t'  •)  ) 

J. 591)38 

-0.92653) 

0.55*565 

1.CC8686 

0.984415 

2.30/145 

J.96’)t)-JJ 

J.  465  753 

-3.0296C2 

”0.596464 

~ T7C04468  

0.5420 32 

2.3218*3 

1.096203 

1. JJC9O0 

- ). 0202.68 

1.C0J9C0 

1.0)034 

l.OCOJJO 

2.336396 

MINE 

in 


PI  TO 

U. 392UJ6 
U. 436770 
0.496236 
0.556329 
0.622756 
0.669720 
3.  759409 

- 0. 831279 
0.90J932 
0.976610 
1.048937 
1.120709 
1. 191607 

- 1.262153 
1.331645 
1.400391 
1.468212 
1.535129 
1.601123 
1.666174 
1.730269 
1.793396 

—1.655554 

1.916735 

1.976944 

2.036166 

~2.09447B“ 

2.151031 

2.208266 

2.263816 

2.316503 

~ ZV372)66“ 
2.42  5444 
2.471779 

“2.529417* 
2.560409 
2. 630806 
2.660668 
2.73UC47 
2.779005 
2.62  7601 
2. 675696 
2.923959 
2.971U47 
3.019625 
3.067357 
3.115106 
3.162937 
3.213911 
3.259093 
3.307544 
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H 

k 

M 

00 

n 

nsh 

P 

q 

r 

R 

s 

St 

t 

T 

* 

u 

* 

Uo= 

U 


B 


SYMBOLS 

viscosity  law  constant,  c*  = 198. 6°R 

* * *2 

skin  friction  coefficient,  2t  / ( u^  ) 

specific  heat  of  constant  pressure 

* * 2 

nondimensional  total  enthalpy,  H /u^ 
thermal  conductivity 
free  stream  Mach  number 

coordinate  measured  normal  to  the  body,  nondimen- 
sionalized  by  the  body  nose  radius 

shock  stand  off  distance  normal  to  the  body  surface 

* * 2 * 

nondimens ional  pressure,  p /(p^  ) 

* * *3 

nondimensional  heat  transfer,  q /( p ^ ) 

nondimensional  axisymmetric  radius 
defined  as  yB  + W05* 

nondimensional  surface  distance  coordinate 
Stanton  number,  St  = qw/(HQ-Hw) 

time,  and  also  used  for  normalized  temperature,  T/T^ 

* *2  *> 

nondimensional  temperature,  T = T /(uB  /Cp) 
free  stream  temperature 

nondimensional  velocity  component  tangent  to  the  body 
surface,  u*/u*  . 

free  stream  velocity 

nondimensional  component  of  velocity  aft  and  tangent 
to  the  shock  interface 

nondimensional  velocity  component  normal  to  the  body 
surface,  v*/u* 

nondimensional  component  of  velocity  aft  and  normal 
to  shock  interface 

axial  distance  for  body  surface  measured  from 
stagnation  point 
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AEDC-TR-79-25 


x . defined  as  x_  - n sin<> 
sh  B s 

yB  normal  distance  for  body  surface  measured  from  axis 
a shock  angle,  see  Figure  1 

6 angle  defined  in  Figure  1 

y ratio  of  specific  heats 

* *2  * * * * 1/2 

e perturbation  parameter,  £ = [y  (uM  /C  ) / P^u^a  ] 

ic  nondimensional  surface  curvature 

* * *2  * 

U nondimensional  coefficient  of  viscosity,  y = U /y  (ua  /C  ) 

* * 

p nondimensional  density,  p = p /p^ 

free  stream  density 

* * *2 

r nondimensional  shear  stress,  t /(p^  } 

body  angle  defined  in  Figure  1 
a Prandtl  number,  a - yC  /K 

E,  nondimensional  surface  distance  coordinate,  s 

g normalized  normal  coordinate,  n/n  ^ 

Subscripts 
1 wall  value 

0 stagnation  conditions 

n used  for  normal  derivatives 

N shock  value 

s used  for  longitudinal  derivatives 

sh  conditions  immediately  behind  the  shock  wave 

°°  free  stream  conditions 

w wall  value 

E,  used  for  longitudinal  derivatives 

ri  used  for  normal  derivatives 
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Superscripts 

• longitudinal  derivative,  d/ds 

physical  quantities  normalised  by  their  shock  values 

* dimensional  quantities,  also  used  for  first  sweep 
of  ADI  numerical  scheme 

j 0 for  plane  flow  and  1 for  axisymmetric  flow 
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